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This  report  describes  a  small  follow-up  effort  to  a  previous  study  of  dynamic  tensile 
failure  of  concrete  in  which  our  long-range  objective  is  to  understand  and  quantify 
the  micromechanics  of  dynamic  tensile  failure  of  concrete.  In  the  previous  study,  we 
developed  an  experimental  technique  to  apply  dynamic  tension  to  5-cm-diameter  x  76-cm- 
long  concrete  rods  at  a  strain  rate  of  about  10/s,  and  we  performed  posttest  computa¬ 
tions  with  a  simple  one-^dimensional  strainrsoftening  model  to  interprete  an  initial 
set  of  experiments.  In  the  current  effort,  our  primary  task  was  to  scrutinize  the 
technique  for  observing  microcracks  in  damaged  specimens  and  to  quantify  the  micro¬ 
cracks  in  some  of  the  specimens  already  tested.  A  second  task  was  to  computationally 
interpret  more  of  the  experimental  results.  Third,  we  made  a  preliminary  step  toward 
computing  the  strength  and  modulus  of  a  material  cell  containing  several  interacting 
cracks . 
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19  ABSTRACT  (Continued) 

Three  methods  for  inspecting  the  concrete  specimens  for  microcracking  were  evaluated. 
The  first  method  was  to  use  a  scanning  electron  microscope  (SEM)  to  view  the  concrete 
surface.  We  found  this  method  to  be  unsatisfactory  because  of  extensive  cracking 
caused  by  evacuating  the  specimen.  The  second  method  was  to  replicate  the  specimen 
surface  with  acetylcellulose  replicating  film  and  to  view  the  film  with  the  SEM. 

This  method  introduces  uncertainties  in  identifying  cracks.  The  third  method  was 
to  view  a  polished  concrete  surface  with  an  optical  microscope  at  a  magnification 
of  lOOx.  This  is  the  preferred  method  for  observing  microcrack  damage  produced  in 
the  dynamic  tension  tests.  ^Microcracks  as  small  as  2  Um  wide  and  100  pm  long  can 
be  seen.  Better  resolution  could  probably  be  attained  with  more  highly  polished 
specimens.  In  the  three  5-cm-long  specimens  inspected,  we  saw  that  microcracks  pass 
through  aggregates  and  around  aggregates,  some  appear  to  be  blunted  by  aggregates, 
and  some  terminate  in  the  mortar.  Nearly  all  of  the  damage  was  found  within  a  dis¬ 
tance  of  J  cm  from  the  primary  fracture.-  Inspection  of  additional  specimens  that 
were  tested  in  dynamic  tension  is  recommended.  •' 

We  used  the  simple  strain-softening  model  to  computationally  interpret  two  additional 
experiments.  Our  approach  was  to  match  the  strain  histories  and  tensile  damage.  The 
experiments  were  performed  on  two  concretes  whose  static  tensile  strengths  are  about 
3.5  MPa.  One  concrete  has  about  twice  the  apparent  dynamic  strength  of  the  other 
concrete.  Sjveral  more  of  the  experiments  performed  previously  should  be  interpreted 
with  the  strain-softening  model  to  extend  our  knowledge  of  how  the  concrete  behaved 
in  these  experiments. 

A  computer  program  was  developed  to  evaluate  Kachanov's  solution  for  stress  intensity 
factors,  crack  face  displacements,  and  effective  moduli  for  an  elastic  material  con¬ 
taining  a  two-dimensional  array  of  cracks.  An  improvement  to  account  for  the  presence 
of  cell  boundaries  is  recommended. 
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SUMMARY 


Background 

This  report  describes  a  small  follow-up  effort  to  a  previous  study  of  dynamic 
tensile  failure  of  concrete  [1,2].  Our  long  range  objective  is  to  understand  and 
quantify  the  micromechanics  of  dynamic  tensile  failure  of  concrete.  In  the  previ¬ 
ous  study,  we  developed  an  experimental  technique  to  apply  dynamic  tension  to 
5-cm-diameter  x  76-cm-long  concrete  rods  at  a  strain  rate  of  about  10/s  [3] 1  and 
performed  posttest  computations  with  a  simple  one-dimensional  strain-softening 
model  to  interpret  an  initial  set  of  experiments  [4]1.  The  strain-softening  compu¬ 
tations  and  a  preliminary  posttest  microscopic  inspection  of  one  of  the  specimens 
indicated  that  this  experiment  produces  distributed  tensile  cracking  damage  that 
can  be  quantified  and  related  to  the  load  history.  We  then  performed  a  larger 
set  of  dynamic  tension  experiments  intended  to  provide  a  range  of  damage  levels 
[2].  The  results  of  these  experiments  and  the  damaged  specimens  are  available  for 
further  analysis. 

In  the  current  effort,  our  objective  was  to  validate  the  experimental  and  ana¬ 
lytical  approach  used  in  the  previous  study.  Our  primary  task  was  to  scrutinize 
the  technique  for  observing  microcracks  in  damaged  specimens  and  to  quantify 
the  microcracks  in  some  of  the  specimens  already  tested.  A  second  task  was  to 
use  a  strain-softening  model  to  computationally  interpret  some  of  the  previous 
experiments.  Finally,  we  made  a  preliminary  step  toward  computing  the  strength 
and  modulus  of  a  computational  cell  containing  several  interacting  cracks.  The 
results  and  conclusions  of  these  tasks  are  summarized  below. 

Previous  Work 

An  experimental  method  was  developed  to  study  the  tensile  failure  of  brittle 
geologic  materials  at  strain  rates  of  approximately  10  to  20/s  [3 j.  In  these  experi- 

1  Copies  of  these  articles  are  attached  as  Appendices  A  and  B. 
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ments,  a  cylindrical  rod  specimen  is  first  loaded  in  static  triaxial  compression,  then 
the  axial  pressure  is  released  from  each  end  simultaneously  and  very  rapidly.  The 
resulting  rarefaction  waves  interact  in  the  center  of  the  rod  to  produce  a  dynamic 
tensile  stress  equal  in  magnitude  to  the  original  static  compression.  The  pressure 
acting  on  the  radial  surface  is  approximately  constant  during  the  experiment.  As 
an  application  of  this  method,  several  experiments  were  performed  on  concrete. 
Transient  measurements  were  made  of  the  axial  load  at  each  end,  the  confining 
pressure,  and  the  axial  and  circumferential  surface  strains  at  several  locations 
along  the  length  of  the  rod. 


Usually  a  single  fracture  occurred  near  the  midpoint  of  the  rod.  In  some 
experiments  multiple  fractures  occurred.  If  we  assume  the  peak  observed  strains 
in  these  experiments  to  be  elastic,  we  estimate  the  unconfined  tensile  strength 
of  the  concrete  at  a  strain  rate  of  10  to  20/s  to  be,  on  average,  approximately 
40%  higher  than  the  static  splitting  tensile  strength.  At  the  same  strain  rate, 
the  tensile  strength  with  10  MPa  confining  pressure  averaged  approximately  100% 
higher  than  the  static  splitting  tensile  strength  and  40%  higher  than  the  unconfined 
tensile  strength  at  10  to  20/s.  Nonlinear  analyses  indicate  that  these  estimates 
are  reasonable  but  that,  in  general,  the  assumption  of  elastic  response  is  not  valid. 
Matching  the  measured  strain  histories  with  calculations  requires  that  the  rod  be 
modeled  inelastically. 


A  one-dimensional  strain-softening  model  was  used  in  wave-propagation  calcu¬ 
lations  to  interpret  the  results  of  the  dynamic  tension  experiments  on  concrete  rods 
[4],  The  model  is  based  on  the  assumption  that  the  stress-strain  relation  is  not 
a  property  of  a  material  point  (as  in  continuum  theory)  but  an  average  property 
of  a  finite  volume  of  material  containing  a  developing  crack  or  failure  plane.  The 
stress-strain  relation  thus  has  associated  with  it  a  finite  dimension,  namely  the 
average  crack  separation  distance.  We  used  this  model  to  simulate  two  dynamic 
unconfined  tension  experiments  and,  by  trial-and-error,  obtained  good  agreement 
with  the  measured  axial  strain  histories  in  both  cases. 


In  addition  to  providing  an  estimate  of  the  dynamic  tensile  properties  of  the 
concrete,  these  calculations  suggest  that  tensile  damage  in  the  concrete  was  dis¬ 
tributed  over  several  centimeters.  Finally,  the  calculations  suggest  that  the  strain 
history  measured  a  few  centimeters  from  the  location  of  fracture  is  primarily  a 
function  of  inelastic  wave  propagation  from  the  fracture  location  to  the  strain 
gage  (through  a  region  of  distributed  tensile  damage)  and  is  less  dependent  on  the 
behavior  of  the  material  right  at  the  fracture. 

Microscopic  Inspection  of  Damaged  Specimens 

Three  methods  for  inspecting  the  concrete  specimens  for  microcracking  were 
evaluated.  The  first  method  was  to  use  a  scanning  electron  microscope  (SEM) 
to  view  the  concrete  surface.  We  found  this  method  to  be  unsatisfactory  because 
of  extensive  cracking  caused  by  evacuating  the  specimen.  The  second  method 
was  to  replicate  the  specimen  surface  with  acetylcellulose  replicating  film  and 
view  the  film  with  the  SEM.  This  method  introduces  uncertainties  in  identifying 
cracks.  The  third  method  was  to  view  a  polished  concrete  surface  with  an  optical 
microscope  at  a  magnification  of  lOOx.  This  is  the  preferred  method  for  observing 
microcrack  damage  produced  in  the  dynamic  tension  tests.  Microcracks  as  small 
as  2  m  wide  and  100  (j. m  long  can  be  seen.  Better  resolution  could  probably  be 
attained  with  more  highly  polished  specimens. 

In  the  three  5-cm-long  specimens  inspected,  we  saw  that  microcracks  pass 
through  aggregates  and  around  aggregates,  some  appear  to  be  blunted  by  aggre¬ 
gates,  and  some  terminate  in  the  mortar.  Nearly  all  of  the  damage  was  found 
within  a  distance  of  3  cm  from  the  primary  fracture.  We  recommend  inspection 
of  additional  specimens  that  were  tested  in  dynamic  tension. 

Strain-Softening  Computations 


site  of  localization  (fracture)  and  that  the  fundamental  property  of  the  material  is 
the  relation  between  stress  and  fracture  volume  per  unit  area  (average  crack  open¬ 
ing).  Within  this  framework,  the  material  cell  dimension  represents  the  spatial 
freqency  of  localization  sites  in  an  inhomogeneous  material.  In  the  previous  study, 
a  0.635-cm  material  cell  size  was  found  to  give  the  best  agreement  with  measured 
strains  and  observed  fractures  in  the  four  experiments  simulated. 

We  used  the  simple  strain-softening  model  to  match  the  strain  histories  and 
tensile  damage  in  dynamic  tension  tests  on  two  concretes  whose  static  tensile 
strengths  are  about  3.5  MPa.  The  same  0.635-cm  material  cell  size  was  chosen, 
and  the  relation  between  stress  and  fracture  volume  per  unit  area  was  adjusted  to 
match  the  experimental  results.  The  apparent  dynamic  strength  of  the  concrete 
used  in  Tests  101  to  106  is  about  twice  as  high  as  the  dynamic  strength  of  the 
concrete  used  in  Tests  41  to  46. 

Several  more  of  the  experiments  performed  previously  should  be  interpreted 
with  the  strain-softening  mode]  to  extend  our  knowledge  of  how  the  concrete  be¬ 
haved  in  these  experiments.  We  also  recommend  additional  study  of  the  apparent 
natural  material  cell  size  for  a  better  understanding  of  its  source  and  meaning. 

Properties  of  a  Multiply  Cracked  Material  Cell 

We  developed  a  computer  program  to  evaluate  a  solution  derived  by  Kachanov 
[7]  for  determining  stress  intensity  factors,  crack  face  displacements,  and  effective 
moduli  for  elastic  material  containing  a  two-dimensional  array  of  cracks.  The 
actual  locations  of  the  cracks  in  the  body,  and  their  interactions  with  each  other 
are  treated.  The  Kachanov  solution  provides  stress  intensity  factors  that  are  within 
10%  of  exact  analytical  solutions  for  cases  with  crack  tip  separations  greater  than 
1%  of  the  crack  length.  Hence,  this  analysis  is  useful  for  studying  crack  interactions 
when  the  cracks  are  fairly  near  each  other;  that  is,  when  the  cracks  are  approaching 
coalescence  and  fragmentation. 


4 


The  limitation  of  the  Kachanov  approach  is  in  the  determination  of  the  effec¬ 
tive  moduli  for  a  finite  cell  of  material,  because  this  approach  does  not  consider 
boundary  conditions  except  those  at  infinity.  A  possible  improvement  to  account 
for  the  presence  of  cell  boundaries  would  be  to  lay  out  a  repeating  set  of  crack 
arrays  so  that  the  computational  cell  being  considered  is  bounded  by  the  symme¬ 
try  planes  between  the  arrays  of  cracks.  We  could  then  construct  the  compliance 


matrix  for  the  central  cell. 
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MICROSCOPIC  INSPECTION  OF  DAMAGED  SPECIMENS 


Specimen  Preparation 

Segments  from  some  of  the  specimens  tested  in  the  previous  project  were  pre¬ 
pared  as  illustrated  schematically  in  Figure  1.  The  5-cm-long  segments  are  cut 
from  the  rod  with  a  diamond-impregnated  saw  blade.  Then  they  are  potted  in  a 
cylinder  of  cement  paste  and  sectioned  along  a  diametric  plane.  The  face  of  the 
sectioned  segment  is  then  polished  with  600-grit  sandpaper. 

Because  most  of  the  microcracking  is  expected  to  have  occurred  within  a  few 
centimeters  of  the  primary  fracture,  we  prepared  specimens  taken  from  this  region 
in  the  rods  used  in  Tests  41  and  103.  Segments  taken  from  the  top  side  of  the 
primary  fracture  are  labeled  ’Tl’  and  ’T2’  and  those  from  the  bottom  side  of  the 
fracture  are  labeled  ’Bl’  and  ’B2’,  with  the  numbers  indicating  the  order  of  the 
segments  from  the  primary  fracture.  The  top  and  bottom  directions  refer  to  the 
orientation  of  the  rod  during  fabrication  and  testing.  The  coordinate  system  used 
for  strain  gage  placement  is  centered  at  the  midpoint  of  the  rod,  with  the  positive 
axis  directed  toward  the  bottom  end  of  the  rod. 

Inspection  Methods 

Three  methods  for  inspecting  the  concrete  specimens  for  microcracking  were 
evaluated.  The  first  method,  using  a  scanning  electron  microscope  (SEM)  to  view 
the  concrete  surface,  was  found  to  be  unsatisfactory  because  of  extensive  cracking 
caused  by  evacuating  the  specimen.  The  second  method,  replicating  the  specimen 
surface  with  acetylcellulose  replicating  film  and  viewing  the  film  with  the  SEM, 
introduces  uncertainties  in  identifying  cracks.  The  third  method,  viewing  the 
concrete  surface  with  an  optical  microscope,  is  the  most  straightforward  and  the 
most  reliable  of  the  three,  but  it  is  limited  to  a  lower  level  of  magnification  than 
methods  using  the  SEM. 

We  first  attempted  to  inspect  the  concrete  for  microcracking  damage  by  view- 
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Figure  1 .  Preparation  of  concrete  specimen  for  microscopic  inspection. 


ing  the  polished  concrete  surface  with  the  SEM.  Unfortunately,  this  method  intro¬ 
duces  extensive  damage  to  the  concrete.  Figure  2  shows  a  pair  of  SEM  photographs 
of  a  virgin  concrete  specimen.  The  network  of  microcracks  results  from  evacuating 
the  specimen  in  preparation  for  coating  it  with  a  conductive  layer  (a  vacuum  is 
also  necessary  for  viewing  with  the  SEM).  We  proved  that  the  cracks  are  caused 
by  the  vacuum  by  inspecting  replicas  of  a  specimen  before  and  after  evacuation 
(not  shown) .  Thus,  even  though  other  researchers  have  utilized  this  or  a  similar 
method  for  microscopic  inspections  of  concrete  [5,6],  we  found  it  unsatisfactory 
for  our  purposes. 

The  second  method  was  to  use  acetylcellulose  replicating  film.  When  applied 
with  a  solvent  to  the  polished  surface  of  the  specimen,  this  film  replicates  the 
surface  by  flowing  into  the  cracks  and  voids.  The  film  is  then  dried,  coated  with  a 
thin  layer  of  gold,  and  viewed  with  an  SEM.  In  the  previous  project,  this  method 
was  used  to  chart  the  microcracks  in  cue  specimen  and  produced  seemingly  good 
results.  However,  upon  further  scrutiny  u'  the  current  study,  we  concluded  that 
there  is  significant  uncertainty  in  identifying  microcracks  with  the  replicating  film. 

Figure  3  compares  an  SEM  photograph  of  an  acetylcellulose  replica  with  an 
optical  photograph  of  the  concrete  specimen.  Note  that  because  the  replica  is 
made  face  down,  the  photographs  are  reversed  images  of  each  other.  The  fine 
white  lines  on  the  replica  were  identified  previously  as  microcracks  from  their 
cracklike  appearance.  The  concrete  specimen  does  not  appear  to  have  cracks  at 
these  locations.  At  a  magnification  of  500x,  a  roughness  indicative  of  surface 
flaking  can  be  seen  with  the  optical  microscope  at  the  locations  of  the  apparent 
cracks,  but  a  true  crack  is  not  visible.  Thus,  it  appears  that  the  replica  accentuates 
these  surface  features  so  that  they  resemble  cracks.  This  can  give  misleading 
information,  and  we  conclude  that  this  method  of  inspection  is  too  unreliable  for 
the  present  purpose. 

The  preferred  method  is  simply  to  view  the  concrete  surface  directly  through 
an  optical  microscope.  Figure  4  shows  a  photograph  of  a  microcrack  at  a  different 
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(a)  Microcrack  network 
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location  in  the  specimen  shown  in  Figure  3.  This  crack  is  about  4  mm  in  length; 
only  a  portion  of  its  total  length  is  shown  in  the  photograph.  Admittedly,  the 
optical  photographs  of  the  concrete  are  not  as  sharp  as  the  SEM  photographs  of 
replicas,  but  when  viewed  through  the  microscope  at  lOOx,  the  concrete  surface 
can  be  seen  with  sufficient  detail  to  identify  without  doubt  cracks  as  small  as  2  (im 
wide  and  100  long.  Larger  cracks  can  be  photographed  clearly  at  lOOx.  When 
the  concrete  surface  is  viewed  at  a  magnification  of  500x  the  surface  roughness  is 
greater  than  the  depth  of  field,  but  by  adjusting  the  focus  while  scanning  a  small 
area  the  viewer  can  garner  additional  detail.  This  method  could  be  improved  if 
the  surface  of  the  specimen  could  be  more  highly  polished  and  if  the  optics  could 
provide  a  greater  depth  of  field. 

Results 

Specimens  from  two  rods  tested  in  the  previous  study  and  a  specimen  from 
a  virgin  rod  were  prepared  and  inspected  for  microcracks.  The  previously  tested 
rods  were  made  of  different  concretes  and  were  preloaded  to  different  levels.  A 
range  of  microcrack  damage  was  observed  in  these  specimens.  No  microcracks 
were  found  in  the  virgin  specimen. 

The  first  specimen  inspected  was  41-Bl.  The  static  splitting  strength  of  this 
concrete  was  about  3.4  MPa.  In  Test  41  the  static  preload  was  a  uniaxial  stress  of 
10.8  MPa,  and  the  rod  fractured  at  -5.0  cm  from  the  midpoint.  This  is  the  same  rod 
that  we  inspected  previously  by  the  replica  method.  However,  in  this  investigation 
we  inspected  the  other  half  of  the  same  segment.  Figure  4  shows  a  photograph  of 
this  specimen,  and  Figure  5  shows  a  map  of  the  full  set  of  microcracks  found.  The 
map  shows  that  there  is  a  concentration  of  damage  about  2.5  cm  from  the  primary 
fracture.  This  is  consistent  with  the  predictions  of  strain-softening  calculations, 
as  will  be  shown  below.  The  long  crack  running  in  the  axial  direction  lies  along 
the  boundary  between  mortar  and  a  long,  slender  aggregate. 

The  next  rod  inspected  was  from  Test  103.  The  static  splitting  strength  of 


this  concrete  was  about  3.6  MPa.  In  Test  103  the  static  preload  was  17.3  MPa; 
the  primary  fracture  occured  at  +2.9  cm  from  the  midpoint.  Figure  6  shows  pho¬ 
tographs  taken  of  specimens  103-Tl  and  103-Bl,  and  the  microcrack  maps  of  these 
specimens  are  shown  in  Figures  7  and  8.  In  103-Tl  the  damage  is  concentrated 
about  2  cm  from  the  primary  fracture;  in  103-Bl  there  is  a  large  crack  about  3 
cm  from  the  primary  fracture.  Again,  strain-softening  calculations  predicted  a 
concentration  of  damage  a  few  centimeters  from  the  primary  fracture. 

Conclusions  and  Recommendations 


Viewing  a  polished  concrete  surface  with  an  optical  microscope  at  a  magni¬ 
fication  of  lOOx  is  a  suitable  method  for  observing  microcrack  damage  produced 
in  the  dynamic  tension  tests.  Microcracks  as  small  as  2  /xm  wide  and  100  /xm 
long  can  be  seen.  Better  resolution  could  probably  be  obtained  with  more  highly 
polished  specimens. 

In  the  three  5-cm-long  specimens  inspected  we  saw  that  microcracks  pass 
through  aggregates  and  around  aggregates;  some  appear  to  be  blunted  by  ag¬ 
gregates;  and  some  terminate  in  the  mortar.  Nearly  all  of  the  damage  was  found 
within  3  cm  of  the  primary  fracture.  We  recommend  inspection  of  additional 
specimens  that  were  tested  in  dynamic  tension  so  that  we  can  relate  quantified 
observations  of  damage  to  known  loading  conditions  for  a  range  of  loads. 
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microscope  photographs  of  polished  concrete  specimens 


COMPUTATIONAL  INTERPRETATION  OF  EXPERIMENTS 


The  purpose  of  the  calculations  performed  in  this  effort  is  to  interpret  the  re¬ 
sults  of  the  experiments.  In  addition  to  the  initial  conditions  and  the  boundary 
conditions,  the  data  available  from  an  experiment  include  the  location  of  fracture 
and  the  strain  histories  at  a  few  locations  on  both  sides  of  the  fracture.  These  data 
do  not,  in  themselves,  provide  a  direct  measure  of  a  material  property.  Our  ap¬ 
proach  to  interpreting  the  results  is  to  make  an  estimate  of  the  material  behavior, 
namely  the  stress-strain  path  followed  in  the  experiment,  and  to  computation¬ 
ally  simulate  the  experiment.  By  adjusting  the  estimate  of  the  stress-strain  path, 
we  can  satisfactorily  match  the  data  from  the  experiment.  Then  the  assumed 
stress-strain  path  is  a  reasonable  estimate  of  the  actual  material  response  in  that 
experiment.  Four  experiments  were  interpreted  in  this  manner  previously.  Two 
more  experiments  were  interpreted  in  the  present  effort. 

Test  41 

Tests  41,  42,  and  43  were  performed  with  a  static  uniaxial  preload  of  about  10.5 
MPa  on  rods  made  of  the  same  concrete.  In  all  three  tests  the  primary  fracture 
was  within  1  cm  of  the  midpoint  of  the  rod.  Tests  42  and  43  were  instrumented 
with  strain  gages,  and  the  results  were  matched  previously  with  computations  with 
a  strain-softening  material  description  [4]  (Appendix  B).  Because  posttest  static 
splitting  tests  were  performed  on  the  rods  from  Tests  42  and  43,  specimens  from 
these  rods  are  not  available  for  inspection  for  microcracks.  Test  41  was  performed 
with  no  instrumentation,  and  it  was  not  previously  simulated.  In  the  current 
effort,  one  5-cm-long  segment  from  this  rod  (41-Bl)  was  inspected  for  microcracks 
(Figures  4  and  5).  We  also  computed  the  response  of  Test  41  using  the  same 
strain-softening  parameters  used  to  simulate  Tests  42  and  43.  We  then  compared 
the  distribution  of  peak  fracture  volume  per  unit  area  (<5)  with  the  map  of  the 
microcracks  observed  in  specimen  41-Bl. 

Figure  9  shows  the  material  description  used  to  computationally  simulate  Test 
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41.  The  a  —  <5  relation,  the  0.635-cm  cell  size,  and  the  a  —  e  relation  are  the  same  as 
used  for  Tests  42  and  43.  The  location  of  the  weak  cell  corresponds  to  the  observed 
location  of  the  primary  fracture  in  the  rod  (+0.7  cm).  In  the  computation,  the 
weak  cell  failed  completely  and  several  cells  on  both  sides  of  the  weak  cell  softened 
without  breaking.  The  computed  peaks  of  the  fracture  volume  per  unit  area  in 
each  cell  corresponding  to  specimen  41-Bl  are  compared  with  the  microcrack  map 
of  that  specimen  in  Figure  10.  The  location  of  the  computed  maximum  tensile 
damage  agrees  with  the  location  of  the  observed  concentration  of  microcracks  to 
within  about  1  cell  width. 


Test  103 


Tests  101  to  106  were  all  performed  on  rods  made  of  a  concrete  similar  to  but 
not  the  same  as  that  used  in  the  earlier  tests  [2].  The  static  uniaxial  preload  in 
Test  103  was  17.3  MPa.  Strain  gages  were  used  at  ±2.54  cm,  ±7.62  cm,  and  ±15.2 
cm  from  the  midpoint  of  the  rod.  Fracture  occurred  at  +2.86  cm. 


The  strain  histories  from  this  test  were  simulated  computationally  using  the 
assumed  material  description  shown  in  Figure  11.  The  strength  is  much  higher 
than  that  used  in  the  simulations  of  Tests  41,  42,  and  43,  but  the  material  cell 
size  and  the  critical  fracture  volume  per  unit  area  are  about  the  same.  Two  types 
of  unloading  were  used.  The  one  with  stress  and  strain  returning  to  zero  implies 
that  all  cracks  close  completely  during  unloading;  the  one  with  the  unloading 
slope  equal  to  the  elastic  loading  slope  implies  that  the  cracks  remain  open  during 
unloading. 


Figure  12  shows  the  comparisons  between  the  measured  axial  strains  and  those 
computed  with  the  assumed  material  properties.  Some  of  the  strain  records  were 
matched  better  with  unloading  to  the  origin;  others  were  matched  better  with 
elastic  unloading.  This  apparent  inhomogeneity  was  also  present  in  earlier  results 
[1.4]. 


COMPUTED  FRACTURE  VOLUME  PER  UNIT  AREA  (pm) 


Figure  1 0.  Comparison  of  computed  and  observed  damage  in  specimen  41  -Bi . 


Figures  13  and  14  compare  the  computed  peaks  of  fracture  volume  per  unit 
area  computed  in  each  cell  corresponding  to  specimens  103-B1  and  103-Tl  with  the 
microcrack  maps  of  these  two  specimens.  The  computed  concentrations  of  damage 
at  about  ±2.5  cm  from  the  primary  fracture  are  consistent  with  the  observations. 
However,  the  predicted  damage  at  about  ±4.5  cm  from  the  primary  fracture  was 
not  seen  in  the  microscopic  inspection.  Perhaps  we  will  find  damage  in  the  adjacent 
specimens  (103-B2  and  103-T2). 

Conclusions  and  Recommendations 

The  strain-softening  model  we  use  to  interpret  the  dynamic  tension  experi¬ 
ments  on  concrete  is  based  on  the  concept  that  a  material  cell  contains  a  single 
site  cf  localization  (fracture)  and  that  the  fundamental  property  of  the  material  is 
the  relation  between  stress  and  fracture  volume  per  unit  area  (average  crack  open¬ 
ing).  Within  this  framework,  the  material  cell  dimension  represents  the  spatial 
freqency  of  localization  sites  in  an  inhomogeneous  material.  In  the  previous  study, 
a  0.635-cm  material  cell  size  was  found  to  give  the  best  agreement  with  measured 
strains  and  observed  fractures  in  the  four  experiments  simulated. 

We  used  the  simple  strain-softening  model  to  match  the  strain  histories  and  ten¬ 
sile  damage  in  dynamic  tension  tests  on  two  concretes  with  static  tensile  strengths 
of  about  3.5  MPa.  We  chose  the  same  material  cell  size  used  previously  (0.635- 
cm)  and  adjusted  the  relation  between  stress  and  fracture  volume  per  unit  area  to 
match  the  experimental  results.  The  apparent  dynamic  strength  of  the  concrete 
used  in  Tests  101  to  106  is  about  twice  as  high  as  the  dynamic  strength  of  the 
concrete  used  in  Tests  41  to  46. 

Several  more  of  the  experiments  performed  previously  should  be  interpreted 
with  the  strain-softening  model  to  extend  our  knowledge  of  how  the  concrete  be¬ 
haved  in  these  experiments.  We  also  recommend  additional  study  of  the  apparent 
natural  material  cell  size  for  a  better  understanding  of  its  source  and  meaning. 


(a)  Microcrack  map. 


(b)  Computed  damage. 


DISTANCE  FROM  PRIMARY  FRACTURE  (cm) 


Figure  13.  Comparison  of  computed  and  observed  damage  in  specimen  103-B1. 
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COMPUTATION  OF  EFFECTIVE  MODULI 
FOR  A  CRACKED  MATERIAL 


A  computational  model  for  representing  the  fracturing  processes  in  concrete 
and  the  mechanical  response  of  concrete  to  loading  must  contain  descriptions  for 
the  initiation,  growth,  and  coalescence  processes  for  cracks.  In  addition,  it  is  nec¬ 
essary  to  know  the  effective  stiffness  of  the  partially  cracked  material  at  all  stages 
of  damage.  The  current  study  of  effective  moduli  is  an  effort  toward  evaluating 
these  stiffnesses. 

Solution  for  Stress  Intensity  Factors  and  Crack  Openings 

Kachanov  [7]  has  constructed  an  approximate  analysis  for  computing  stress  in¬ 
tensity  values  for  arrays  of  cracks  of  arbitrary  locations,  lengths,  and  orientations 
in  an  elastic  solid.  The  method  is  based  on  classical  solutions  for  the  stress  states 
around  a  two-dimensional  flat  crack  in  an  infinite  elastic  body  under  external  trac¬ 
tions.  The  procedure  provides  for  a  simultaneous  solution  for  the  stress  intensity 
values  (FT/  and  Ku )  at  each  end  of  each  crack  in  the  array.  With  these  K  values, 
the  crack  opening  shape  is  determined.  From  the  crack  shape,  the  average  crack 
opening  strain  in  the  body  can  be  found. 

Kachanov  considers  an  array  of  two-dimensional  cracks  embedded  in  an  elastic 
material  under  a  remote  loading  (ct00).  The  crack  lengths  are  given  by  L,-  and  the 
orientations  by  the  unit  normals  n<. 

The  stress  acting  along  each  crack  is  a  superposition  of  the  external  stress  field 
on  the  stress  field  caused  by  the  opening  of  each  of  the  other  cracks.  To  proceed, 
Kachanov  considers  just  two  cracks  and  their  interactions.  He  evaluates  first  the 
interaction  stresses  cr,"(£)  and  a£(£):  the  stresses  at  position  £  along  the  jth  crack 
generated  by  uniform  normal  and  shear  tractions  of  unit  intensity  along  the  ith 
crack.  These  stress  quantities  are  obtained  from  a  standard  elastic  analysis  of  the 
stresses  in  an  infinite  body  resulting  from  normal  and  shearing  stresses  applied 


to  the  faces  of  an  embedded  crack  (here  the  ith  crack).  To  obtain  the  actual 
stresses  along  the  jth  crack,  we  multiply  cr,”(£)  and  cr,?(£)  by  the  average  of  the 
normal  or  shear  stresses  on  the  ith  crack:  <  Pi  >  and  <  r,-  >.  This  use  of  the 
average  stresses  on  the  cracks  (instead  of  the  actual  strongly  varying  stresses)  is 
the  essential  ingredient  of  Kachanov’s  method  that  simplifies  the  solution  enough 
to  make  it  practical. 

After  the  interaction  stresses  are  determined,  the  stress  state  along  each  crack 
is  computed.  From  this  stress  state,  the  Kj  and  Kji  stress  intensity  factors  are 
evaluated  for  each  crack. 

Next,  these  stress  solutions  are  used  to  evaluate  the  normal  opening  of  the 
crack  and  its  shearing  displacement.  These  procedures  are  given  in  more  detail 
in  Appendix  C.  A  computer  program,  KCRACK,  was  written  for  routine  evalua¬ 
tion  of  stress  intensity  factors  and  crack  opening.  The  solution  procedure  is  also 
described  in  Appendix  C  and  the  program  is  listed  there. 

Results  are  given  in  Appendix  C  for  stress  intensity  factors  for  pairs  of  identical 
cracks  along  the  x-axis.  The  K  values  are  determined  as  a  function  of  the  separa¬ 
tion  of  the  cracks  and  compared  with  exact  solutions.  The  error  in  the  K  values 
increases  as  the  separation  distance  between  the  cracks  decreases.  The  error  is  less 
than  1%  until  the  crack  tip  separation  is  less  than  10%  of  the  crack  lengths  [7]. 
At  this  point  the  K  value  is  45%  larger  than  the  K  value  computed  for  an  isolated 
crack.  Therefore,  it  appears  that  Kachanov’s  procedure  is  essential  for  evaluating 
the  cracked  state  of  a  body  if  the  crack  tips  are  near  each  other. 

Computation  of  Effective  Moduli 

When  the  crack  stress  intensity  factors  and  the  crack  openings  have  been  com¬ 
puted,  it  is  possible  to  determine  the  effective  moduli  for  a  cracked  body.  Such 
moduli  are  termed  “effective”  because  they  provide  an  average  of  the  response 
of  the  intact  elastic  material  around  the  cracks  plus  the  response  of  the  cracks 


themselves. 

Before  exploring  the  two-dimensional  problem  of  computing  moduli  using  Kacha 
nov’s  method,  we  introduce  the  determination  of  the  longitudinal  modulus  of  a  rod 
of  elastic  material  with  transverse  cracks.  This  rod  problem  will  serve  to  indicate 
the  essence  of  the  method  and  assumptions  behind  the  more  complex  derivation 
of  Kachanov.  Let  us  examine  the  longitudinal  stiffness  of  a  rod  with  radius  R  and 
length  L  with  N  cracks  all  located  such  that  the  normals  to  their  planes  are  along 
the  axis  of  the  rod.  The  cracks  have  radii  iZj.  Under  a  longitudinal  stress  (cr), 
their  tensile  opening  is  given  by 


where  is  one-half  the  maximum  separation  of  the  crack  faces  and  E  and  u 
are  Young’s  modulus  and  Poisson’s  ratio.  The  crack  faces  form  an  ellipsoid  with 
three  semiaxes,  6,,  iZ,,  and  iZ,.  Therefore,  the  volume  of  the  opening  of  the  crack 
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To  determine  an  effective  elongation  (6^)  for  the  crack,  we  average  this  crack 
volume  over  the  area  (7riZ2)  of  the  rod. 
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Now  consider  the  elongation  that  occurs  to  the  rod  when  a  stress  (a)  is  applied. 
We  assume  that  the  intact  elastic  material  is  all  under  the  same  tensile  stress  (cr) 
as  it  would  be  for  an  uncracked  rod,  so  the  elongation  is 
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In  addition  to  this  elongation,  we  have  the  elongation  given  by  all  the  cracks: 
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The  effective  modulus  (Me)  simply  relates  the  applied  stress  (a)  to  the  average 
strain  induced  in  the  rod. 


Me  = 


(ALm  +  ALc)/L  ,  16(1  —  i/2 

1  l  0 1/ 


where  V,  is  the  volume  of  the  rod:  nR2L.  From  the  equation  for  Me,  we  see 
that  the  stiffness  of  the  cracked  rod  is  reduced  by  the  factor  in  the  denominator 
that  is  a  function  of  the  cube  of  the  crack  radii. 


In  deriving  the  effective  moduli  for  a  cracked  body,  Kachanov  [7,8]  uses  his 
solution  for  the  stress  intensity  factor  and  crack  opening  for  an  array  of  cracks. 
He  computes,  one  at  a  time,  the  compliance  factors  ( Ciju ),  where 


C«;  =  Cijkl°kl 


and  €,,■  and  ou  are  the  average  strain  and  stress  on  the  cracked  body.  The 
compliance  term  CfJ*i  is  obtained  by  applying  the  stress  <tju  to  the  body,  computing 
all  the  strain  and  distortion  quantities  for  the  body  and  cracks,  and  summing 
those  which  contribute  to  e,;.  Ciju  is  then  computed  by  inverting  the  preceding 


equation.  The  corresponding  effective  stiffness  matrix  is  computed  by  inverting 
the  compliance  matrix.  Details  of  this  derivation  are  given  in  Appendix  C. 

With  a  procedure  available  for  computing  the  moduli  for  a  cracked  body,  we 
can  undertake  a  number  of  computations: 

•  Evaluate  moduli  for  some  cases  for  which  there  are  analytical  solutions. 

•  Compute  the  moduli  for  simple  parallel  arrays  of  cracks  and  compare  these 
results  with  the  results  of  simpler  procedures. 

•  Evaluate  the  moduli  for  arrays  of  cracks  and  determine  the  sensitivity  of  the 
moduli  to  closeness  of  the  cracks. 

No  analytical  solutions  were  found  that  matched  the  cases  possible  with  the 
Kachanov  method,  so  the  first  type  of  computation  was  not  performed.  The 
other  types  of  studies  were  undertaken,  and  the  results  are  given  below. 

Parallel  Cracks  in  Line.  The  first  analysis  concerns  a  pair  of  horizontal  cracks 
of  length  1  —  a  with  a  spacing  of  2 a  between  their  nearer  tips.  A  stress  (auv)  was 
applied  to  a  block  containing  these  cracks.  The  results  in  Figure  15  show  the  stress 
intensity  factors  and  compliance  (C22)  for  a  range  of  a  values  and  cross-sectional 
areas  (5).  The  relative  stress  intensities,  Ki/Kjq  (stress  intensity  at  the  near  tip 
normalized  by  the  stress  intensity  for  an  isolated  crack) ,  increase  markedly  during 
this  interval  of  a  values.  The  compliance  factors  increase  more  moderately  over 
the  same  range  of  o.  These  compliances  are  normalized  such  that  the  compliance 
for  a  crack-free  body  is  1.  Here  the  results  for  C22  are  grouped  according  to  the 
area  (5).  The  effect  of  the  area  can  be  easily  accounted  for  as  follows: 


Cs=30  —  — (Cs=2  —  1)  +  1 


Therefore,  the  area  was  held  constant  ( S  =  2)  in  the  rest  of  the  study. 
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Figure  1 5.  Variation  of  compliance  and  normalized  stress  intensity  with  crack 
separation  for  two  cracks  in  line. 

S  =  area. 


Pairs  of  Cracks,  One  Above  the  Other.  The  next  series  of  calculations  was 
made  for  pairs  of  parallel  cracks,  one  above  the  other,  with  a  spacing  (6).  The 
results  shown  in  Figure  16  indicate  that  neither  the  stress  intensity  factor  nor 
the  compliance  factor  (C22)  is  very  sensitive  to  the  crack  spacing  (6).  For  large 
spacings,  the  Kj  value  would  certainly  approach  Ki0,  the  value  for  an  individual 
crack.  Hence,  the  presence  of  a  parallel  crack  tends  to  reduce  the  stress  intensity 
factor.  Similarly,  the  compliance  is  reduced  when  the  cracks  are  near  each  other, 
although  the  effect  is  not  large.  The  compliance  for  the  pair  of  cracks  approaches 
that  for  a  single  crack  (1.318  for  the  present  geometry),  as  expected. 

Four  Crack  Sets.  The  third  case  concerns  a  set  of  four  cracks,  two  horizontal 
and  two  vertical,  extending  directly  out  from  the  origin  (as  shown  in  Figure  17). 
For  large  distances  (a)  between  the  crack  tips,  the  response  of  this  array  is  much 
like  that  shown  in  Figure  15  for  two  horizontal  cracks;  but  as  the  spacing  decreases, 
the  vertical  pair  of  cracks  begins  to  influence  the  interaction.  Although  the  stress 
is  applied  in  the  crvv  direction  only,  horizontal  stresses  act  on  the  vertical  cracks, 
causing  a  large  stress  intensity,  K2  (for  this  configuration  K20,  the  stress  intensity 
for  a  noninteract  g  crack,  is  zero).  As  the  spacing  decreases,  this  stress  intensity 
of  the  vertical  cracks  increases  faster  than  the  stress  intensity  of  the  horizontal 
cracks.  The  compliance  (C22)  for  the  noninteracting  case  is  the  same  as  that 
shown  in  Figure  15.  From  the  symmetry  of  the  problem,  here  Cn  =  C22,  with 
or  without  crack  interaction.  Because  of  the  importance  of  the  vertical  cracks  in 
the  effects  on  the  stress  field,  the  compliance  for  interacting  cracks  increases  faster 
with  decreasing  spacing  in  this  case  than  in  Figure  15. 

Conclusions  and  Recommendations 

From  the  simulations  and  the  attempts  to  examine  cases  with  known  results, 
some  of  the  limitations  of  the  Kachanov  approach  became  clear.  These  limita¬ 
tions  are  important  in  the  determination  of  the  effective  stiffness,  but  not  in  the 
determination  of  stress  intensity  factors  or  crack  face  displacements.  The  anal- 
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Figure  1 6.  Variation  of  compliance  and  normalized  stress  intensity  with  vertical 
spacing  between  two  horizontal  cracks. 

S  -  area. 


ysis  treats  in  detail  the  interaction  between  two  or  more  cracks,  but  it  does  not 
consider  boundary  conditions  except  those  at  infinity.  Hence,  one  cannot  readily 
take  a  finite  region  such  as  a  computational  cell,  lay  out  an  array  of  cracks,  and 
determine  the  effective  moduli  of  this  assemblage.  The  size  of  the  computational 
cell  only  enters  the  analysis  with  the  factor  S,  not  with  an  explicit  location  of  cell 
boundaries. 

A  possible  improvement  to  account  for  the  presence  of  cell  boundaries  would 
include  the  following  steps: 

•  Lay  out  a  set  of  crack  arrays  such  that  they  repeat.  For  example,  the  set  in 
Figure  17  could  be  arranged  with  identical  sets  above,  below,  and  on  both 
sides  of  the  set  of  interest.  Then  the  computational  cell  being  considered  is 
bounded  by  the  symmetry  planes  between  the  sets  of  cracks. 

•  Perform  the  analysis  as  above  to  determine  the  stress  intensity  factors  and 
crack  openings  for  all  cracks. 

•  Construct  the  compliance  matrix  for  only  the  central  cell,  using  only  the 
cracks  present  in  that  cell. 

The  foregoing  procedure  would  provide  a  closer  approximation  to  the  compliance 
of  a  bounded  region.  Yet  this  procedure  also  suffers  from  the  assumption  that 
the  compliance  of  the  matrix  material  around  the  crack  is  under  a  uniform  stress 
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An  experimental  method  was  developed  to  study  the  tensile  failure  of  brittle  geologic  materials  at  strain  rates  of 
approximately  10  to  20/s.  In  these  experiments,  a  cylindrical  rod  specimen  is  first  loaded  in  static  triaxial  compression,  then 
the  axial  pressure  is  released  from  each  end  simultaneously  and  very  rapidly.  The  resulting  rarefaction  waves  interact  in  the 
center  of  the  rod  to  produce  a  dynamic  tensile  stress  equal  in  magnitude  to  the  original  static  compression.  The  pressure  acting 
on  the  radial  surface  is  approximately  constant  during  the  experiment.  As  an  application  of  this  method,  several  experiments 
were  performed  on  concrete  Transient  measurements  were  made  of  the  axial  load  at  each  end.  the  confining  pressure,  and 
axial  and  circumferential  surface  strains  at  several  locations  along  the  length  of  the  rod.  Usually  a  single  fracture  occurred 
near  the  midpoint  of  the  rod.  In  some  experiments  multiple  fractures  occurred.  Assuming  the  peak  observed  strains  in  these 
experiments  to  be  elastic,  the  unconfined  tensile  strength  of  the  concrete  at  a  strain  rate  of  10  to  20/s  was.  on  average, 
approximately  40?  higher  than  the  static  splitting  tensile  strength.  At  the  same  strain  rate,  the  tensile  strength  with  10  MPa 
confining  pressure  averaged  approximately  100?  higher  than  the  static  splitting  tensile  strength  and  40?  higher  than  the 
unconfined  tensile  strength  at  10  to  20/s.  Nonlinear  analyses  indicate  that  these  estimates  are  reasonable,  but  that  in  general 
the  assumption  of  elastic  response  is  not  valid.  To  match  the  measured  strain  histories  with  calculations  requires  that  the  rod 
be  modeled  melastically. 


1.  Introduction 

The  study  of  dynamic  tensile  failure  in  geologic 
materials  and  concrete  is  important  to  many  en¬ 
gineering  applications,  such  as  rapid  excavation, 
in-situ  fracture,  and  impulsive  loading.  It  is  also  of 
fundamental  interest  in  the  field  of  mechanics  of 
materials.  Tensile  failure  in  these  materials  is  pro¬ 
duced  by  the  nucleation,  growth,  and  coalescence 
of  microcracks,  and  the  tensile  strength  is  the 
stress  at  which  this  process  of  accumulating 
damage  becomes  locally  unstable.  Most  geologic 
materials  and  concrete  are  brittle,  that  is,  the 
damage  growth  and  strength  reduction  occur  in  a 
very  short,  but  finite  time.  Details  of  these 
processes  are  important  in  applications  in  which 
the  load  duration  is  comparable  to  the  time  re¬ 
quired  for  tensile  failure.  In  such  problems,  there 
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is  a  strong  interaction  between  the  applied  load 
and  the  growth  of  damage.  Hence,  an  important 
need  in  the  study  of  dynamic  tensile  failure  in 
brittle  geologic  materials  is  to  characterize  the 
failure  process  for  a  wide  range  of  strain  rates. 

Tensile  failure  at  strain  rates  greater  than  10?/s 
has  been  observed  in  concrete  (Gupta  and  Sea¬ 
man,  1979)  and  rock  (Grady  and  Kipp,  1979)  in 
uniaxial  strain  plate  impact  experiments.  Kipp. 
Grady,  and  Chen  (1980)  have  shown  that,  at  these 
very  high  strain  rates,  the  tensile  strength  in¬ 
creases  with  strain  rate  to  the  power  of  } .  Tensile 
failure  at  strain  rates  between  10/s  and  100/s  has 
been  produced  by  reflecting  a  compressive  stress 
pulse  from  the  free  end  of  a  long  rod  specimen. 
Birkimer  (1968)  and  Goldsmith,  Kenner,  and 
Richetts  (1968)  used  this  method  in  experiments 
on  concrete,  in  which  the  compressive  pulse  was 
produced  by  an  impact  of  a  spherical  pellet.  Ab¬ 
bott  and  Cornish  (1965)  and  Felix  (1977)  studied 
ceramics  and  oil  shale,  respectively,  using  explo- 
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sive  charges  to  produce  the  compressive  pulse.  At 
these  intermediate  strain  rales  the  tensile  strength 
dependence  on  strain  rate  is  not  understood,  al¬ 
though  various  emperical  models  have  been  de¬ 
scribed  (Reinhardt,  1985).  Part  of  the  difficulty  is 
that  rod  impact  experiments  are  difficult  to  analyze 
because  loading  conditions  are  not  well  de¬ 
fined.  An  accurate  specification  of  the  loading 
conditions  is  needed  for  a  quantitative  analysis  of 
the  data. 

The  effects  of  loading  path,  including  confining 
stresses,  are  also  important  in  studies  of  dynamic 
tensile  failure  of  geologic  materials.  These  effects 
have  received  minimal  study,  however,  because 
loading  rates  and  loading  paths  cannot  be  varied 
independently  with  most  experimental  methods. 
In  each  of  the  studies  mentioned  above,  only 
uniaxial  strain  loading  or  uniaxial  stress  loading 
was  considered. 

In  this  paper  we  describe  the  development  of 
an  experimental  method  that  addresses  the  issues 
indicated  above  for  brittle  geologic  materials.  The 
experiment  is  based  on  the  concept  of  the  interac¬ 
tion  of  two  rarefaction  waves  to  produce  tensile 
failure  in  a  rod-shaped  specimen.  This  experiment 
meets  the  following  objectives: 

(1)  to  measure  the  tensile  response  at  strain 
rates  between  10/s  and  100/s,  including  the 
strength  reduction  with  the  growth  of  tensile  dam¬ 
age, 

(2)  to  examine  the  effect  of  confining  pressure 
on  tensile  failure  at  these  strain  rates,  and 

(3)  to  accurately  characterize  the  loading  con¬ 
ditions. 

As  an  application  of  the  method,  we  have  cho¬ 
sen  to  study  concrete  because  of  the  current  inter¬ 
est  in  its  dynamic  properties,  and  because  of  its 
similarities  with  geologic  solids.  The  method  is 
also  applicable  to  hard  geologic  solids  that  are 
weak  in  tension.  Several  experiments  were  per¬ 
formed  on  concrete  and  selected  results  are  pre¬ 
sented  here.  Complete  details  of  the  experiments 
and  accompanying  analyses  are  presented  else¬ 
where  (Gran,  1985). 


2.  Experimental  concept 

Approach 

The  concept  underlying  the  experiments  is  to 
use  the  interaction  of  rarefaction  waves  to  produce 
dynamic  tensile  'tresses  in  a  rod-shaped  specimen. 
A  schematic  view  of  the  experimental  technique  is 
shown  in  Fig.  1.  A  cylindrical  rod  is  initially  held 
in  static  compression,  in  both  the  axial  and  radial 
directions.  The  axial  (/>,)  and  radial  (P2)  pres¬ 
sures  are  controlled  separately.  The  radial  pressure 
is  approximately  constant  during  the  experiment. 
The  axial  pressures  at  each  end  of  the  rod  are 
released  simultaneously,  sending  axial  rarefaction 
waves  toward  the  center.  Individually,  these  waves 
bring  the  rod  only  to  zero  axial  stress,  but  when 
they  superpose  at  the  midpoint,  they  bring  the  rod 
to  a  tensile  stress  equal  in  magnitude  to  the  origi¬ 
nal  axial  compression.  Tensile  failure  occurs  near 
the  center  of  the  rod  if  the  resultant  tensile  stress 
exceeds  the  tensile  strength  for  these  conditions. 

Transient  measurements  are  made  of  the  axial 
load  at  each  end,  as  well  as  the  confining  pressure 
and  axial  and  circumferential  surface  strains  at 
several  locations  along  the  length  of  the  rod.  The 
pressure  measurements  define  the  boundary  con¬ 
ditions.  The  strain  measurements  record  the  ef¬ 
fects  of  the  tensile  failure  process  on  the  stress 
waves  propagating  in  the  rod. 

Because  the  specimen  is  initially  loaded  in  static 
axial  compression,  this  method  is  applicable  to 
materials  for  which  the  dynamic  tensile  strength  is 
lower  than  the  static  compressive  strength. 


Specimen  r 
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Fig.  1.  Dynamic  tensile  loading  technique  with  confinemeni 
End  pressures  (F,)  are  released  rapidk  and  simultaneously 
radial  pressure  ( P2 )  is  held. 
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Wave  propagation  preliminaries 

The  rarefaction  waves  propagating  from  the 
ends  of  the  rod  are  not  step  waves  because  the 
applied  pressure  decays  to  zero  in  a  finite  time. 
This  fact  makes  it  more  difficult  to  visualize  the 
stress  distributions  in  the  region  of  wave  interac¬ 
tion.  To  illustrate  and  understand  the  complexities 
introduced  by  the  finite  decay  time,  two  idealized 
examples  of  rod  response  are  considered.  The  first 
example  is  an  elastic  rod  with  no  failure.  The 
second  example  is  a  rod  that  fractures  at  a  single 
point  but  remains  elastic  everywhere  else.  This 
example  can  also  be  visualized  as  two  elastic  half¬ 
rods  connected  together  by  a  weak  bond. 

For  one-dimensional  linear  elastic  response,  the 
axial  stress  or  strain  distribution  along  the  length 
of  the  rod  at  any  time  is  the  solution  of  an  initial 
value  problem  for  the  one-dimensional  wave  equa¬ 
tion.  The  well-known  solution  of  this  problem  is 
the  sum  of  two  waves  traveling  in  opposite  direc¬ 
tions  at  the  same  speed  (Berg  and  McGregor, 
1966).  The  solutions  for  the  following  examples 
were  obtained  using  the  principle  of  superpostion. 

Plots  of  the  axial  stress  and  strain  distributions 
for  several  times  are  shown  in  Fig.  2  for  an  elastic 
rod.  For  simplicity  in  this  example,  a  linear  decay 
in  applied  pressure  is  assumed.  Tensile  stress  is 
taken  to  be  positive.  The  dashed  lines  represent 
the  two  rarefaction  waves  and  the  static  compres¬ 
sion.  The  solid  line  is  the  total  stress  produced  by 
the  superposition  of  the  two  rarefaction  waves  and 
the  preload.  The  distance  the  wave  travels  during 
the  time  required  for  the  pressure  to  drop  to  zero 
is  denoted  by  X.  For  the  assumption  of  one-di¬ 
mensional  response  to  be  appropriate,  X  should 
be  greater  than  a  few  rod  diameters.  In  the  limit¬ 
ing  case  of  instantaneous  pressure  drop,  the  rare¬ 
faction  waves  would  be  step  waves  and  X  would 
be  zero.  The  other  terms  in  the  figure  are  defined 
as  follows:  a0  is  the  magnitude  of  the  peak  axial 
stress,  c0  is  the  corresponding  peak  axial  strain, 
and  X  is  the  position  along  the  length  of  the  rod, 
with  the  origin  at  the  midpoint. 

The  sequence  of  stress  distributions  illustrates 
that  tension  first  occurs  in  a  central  region  of 
width  X,  at  the  time  when  the  rarefaction  waves 
have  overlapped  by  ^X.  As  the  waves  overlap 


Rarefaction  Wave 


jA-4451  -17C 


Fig.  2.  Axial  stress  and  strain  distributions  in  a  one-dimen¬ 
sional  elastic  rod  with  ramp  unloading.  Time  increases  from  (a) 
to  (h). 


further,  the  region  of  tension  broadens.  The  mid¬ 
point  is  the  first  point  to  attain  the  maximum 
tensile  stress,  o0.  However,  intermediate  values  of 
tension  are  attained  simultaneously  in  a  finite 
region.  As  the  tension  increases  to  o0.  the  region 
of  uniform  tension  becomes  narrower.  When  the 
wavefronts  overlap  by  X,  the  maximum  tension  is 
o0  and  exists  only  at  the  midpoint.  Thereafter,  the 
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Fig.  3.  Stress  and  strain  histones  in  a  one-dimcnsiona!  elastic 
rod  with  ramp  unloading. 


region  of  maximum  tension  expands. 

For  this  example,  stress  and  strain  histories  at 
several  points  are  plotted  in  Fig.  3.  This  figure 
shows  that  the  stress  history  is  different  at  every 
point  (except  pairs  of  points  symmetrically  located 
with  respect  to  the  origin).  However,  points  that 
simultaneously  reach  intermediate  values  of  ten¬ 
sion  also  experience  the  same  tensile  stress  history 
(and  tensile  strain  rate)  up  to  that  time.  For 
example,  all  points  between  ±  ^X  reach  a  tension 
of  jo0  simultaneously.  They  also  have  the  same 
stress  history  from  the  time  they  reach  zero  stress 
until  the  time  they  reach  ^o0. 

The  average  tensile  strain  rate — the  peak  tensile 
strain  divided  by  the  rise  time  from  zero  strain  — 
varies  with  position.  The  average  tensile  strain  rate 
for  X  &  is  that  produced  by  a  single  rarefac¬ 


tion  wave  because  the  wavefronts  do  not  overlap 
in  this  region.  For  0  <  X  <  ^X  the  wavefronts  over¬ 
lap,  but  the  average  tensile  strain  rate  varies  with 
distance  from  the  origin  because  the  times  of 
arrival  of  the  two  wavefronts  are  staggered  by 
different  amounts.  The  average  tensile  strain  rate 
at  the  midpoint  is  produced  by  the  simultaneous 
arrival  of  both  rarefaction  waves,  and  is  twice  the 
strain  rate  at  the  more  remote  locations.  Thus,  the 
midpoint  is  the  first  point  to  reach  the  maximum 
stress  and  it  experiences  the  highest  average  tensile 
strain  rate. 

This  example  shows  that  even  without  tensile 
failure,  the  stress  profiles  and  stress  histories  are 
not  simple.  If  o„  exceeds  the  dynamic  tensile 
strength,  the  rod  will  fracture  somewhere  within 
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Fig.  4.  Stress  distributions  in  an  ideally  brittle  rod  with  ramp 
unloading.  Time  increases  from  (a)  to  (f). 
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Fig.  5  Stress  histories  in  an  ideally  brittle  rod  with  ramp 
unloading. 


the  finite  region  that  reaches  the  critical  stress 
simultaneously,  making  subsequent  stress  profiles 
and  stress  histories  even  more  complicated. 

To  illustrate  the  effect  of  fracture,  an  example 
with  fracture  occurring  at  the  midpoint  is  il¬ 
lustrated  in  Figs.  4  and  5.  In  this  example,  the 
fracture  is  assumed  to  be  ideally  brittle  and  to 
occur  instantaneously  when  a  critical  stress  is  re¬ 
ached.  It  is  also  assumed  that  the  rod  remains 
elastic  at  every  point  except  the  midpoint.  The 
critical  stress  for  fracture  is  assumed  to  be  aa0. 
where  0  <  a  <  1.  The  plots  are  drawn  for  a  = 
The  dashed  lines  are  the  static  preload,  the  rare¬ 
faction  waves,  the  reflection  of  one  rarefaction 


wave  from  the  new  free  surface,  and  the  step  wave 
required  to  satisfy  the  new  stress-free  boundary 
condition.  The  solid  line  is  the  total  stress  ob¬ 
tained  by  summation. 

Figure  4(a)  shows  the  stress  distribution  just 
before  fracture.  Figures  4(b)  through  4(f)  show  the 
stress  distributions  after  fracture.  Although  frac¬ 
ture  is  assumed  to  occur  at  the  origin  in  this 
example,  the  figure  illustrates  that,  at  the  time  of 
fracture,  the  critical  stress  exists  throughout  the 
central  region  of  width  A(1  -  a).  When  fracture 
occurs,  the  stress  at  the  midpoint  immediately 
drops  to  zero,  but  the  stress  at  the  other  points  in 
the  central  region  continues  to  increase.  This  is 
because  the  stress  wave  emanating  from  the  frac¬ 
ture  propagates  at  the  same  velocity  as  all  the 
other  waves.  Until  the  effect  of  the  fracture  propa¬ 
gates  to  a  given  point,  the  stress  at  that  point 
continues  to  increase  as  if  fracture  has  not  oc¬ 
curred.  The  portion  of  the  rarefaction  wave  that 
has  already  propagated  past  the  midpoint  by  the 
time  fracture  occurs  has  a  peak  of  ](1  +  a)o0. 
which  exceeds  aa0  for  a  <  1.  That  is.  if  no  other 
fracture  occurs,  the  stress  at  every  point  except  the 
midpoint  will  exceed  the  critical  stress. 

The  stress  histories  at  several  points  in  this 
example  are  plotted  in  Fig.  5.  The  peak  stress  at 
the  midpoint  is  less  than  the  peak  stress  at  every 
other  point.  The  average  :  ensile  strain  rate  at  the 
midpoint  is  equal  to  or  greater  than  the  average 
tensile  strain  rate  at  every  other  point.  The  maxi¬ 
mum  tensile  stress  §(1  +  a)o0  occurs  first  at  ^A. 
The  average  tensile  strain  rate  at  this  point  is  half 
that  at  the  midpoint.  The  tensile  stress  history  is 
different  at  every  point  between  the  origin  and  1  A. 
The  tensile  stress  history  at  all  other  points  is  the 
same  as  at  ^A. 

The  preceding  examples  are  oversimplifications 
of  the  response  of  real  materials  in  this  type  of 
experiment,  but  they  provide  an  insight  into  more 
complex  situations.  The  stress  histories  and  stress 
distributions  in  these  examples  indicate  that  there 
are  a  variety  of  possible  load  histories  leading  to 
failure  and  suggest  that  multiple  fracture  is  likely 
to  occur  in  some  cases.  Experiments  in  which  only 
one  fracture  occurs  are  desirable  because  interpre¬ 
tations  of  them  are  much  more  straightforward. 
The  type  of  load  to  failure  and  the  likelihood  of 
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multiple  fracture  depends  on  the  magnitude  of  the 
preload  in  relation  to  the  tensile  strength  and  the 
rate  dependence  of  the  tensile  properties  of  the 
rod.  In  experiments  with  a  preload  that  is  only 
slightly  higher  than  that  required  to  produce  dy¬ 
namic  failure  during  the  nse  of  the  tensile  stress, 
the  width  of  the  region  that  first  reaches  the 
critical  stress  is  minimized  and  so  is  the  likelihood 
of  multiple  fracture. 


3.  Experimental  method 

Apparatus 

Figure  6  shows  a  photograph  of  the  apparatus 
developed  to  perform  experiments  of  the  type  just 
discussed.  This  apparatus  tests  5.1  cm-diameter, 
76.2  cm-long  rods  at  stresses  up  to  20  MPa.  The 
static  end  pressure  is  removed  in  about  30  ps, 
producing  unloading  strain  rates  in  concrete  of 
about  10/s. 

The  essential  component  of  the  tensile  testing 
apparatus  is  the  unloading  device  at  each  end  of 
the  rod.  Its  design  is  shown  in  Fig.  7.  It  consists  of 
a  bored  aluminum  block  into  which  the  end  of  a 
rod  specimen  and  a  plastic  piston  fit  to  form  a 
chamber  for  oil.  The  oil  is  pressurized  by  means  of 
a  small  orifice  through  the  wall  of  the  bore  block. 
The  rod  and  piston  seal  the  pressurized  oil  in  the 
chamber  with  rubber  O-rings.  The  piston  is  held 
in  place  by  a  thin-walled  steel  support  tube,  which 
presses  against  a  reaction  plate  bolted  to  the  block. 
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Fig.  7.  Unloading  device. 
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The  support  tube  comprises  three  sections,  one  of 
which  is  a  segmented  ring  that  is  explosively  driven 
inward  to  free  the  piston  and  initiate  decompres¬ 
sion  of  the  oil.  The  rod  is  held  in  place  by  an 
identical  unit  at  the  other  end,  with  the  two  blocks 
bolted  together. 
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The  rings  in  the  support  tubes  are  removed 
with  an  estimated  simultaneity  of  less  than  5  ^s 
using  strands  of  sheet  explosive.  The  explosive 
produces  no  damage  to  the  apparatus  except  to 
the  rings,  and  allows  a  smooth  decay  of  pressure 
in  the  oil  chamber.  Occasionally  a  small  compres¬ 
sive  precursor  is  produced  by  the  explosive  pres¬ 
sure. 

In  experiments  without  confinement,  a  Plexi¬ 
glas  tube  is  used  to  space  the  two  unloading 
devices,  as  shown  in  Fig.  6.  In  experiments  with 
confinement,  a  7.6  cm-ID  aluminum  tube  is  used 
to  hold  the  confining  pressure  and  to  space  the 
unloading  devices.  The  confining  pressure  remains 
approximately  constant  during  an  experiment,  but 
it  is  perturbed  slightly  by  the  the  radial  motion  of 
the  rod  produced  by  the  stress  waves. 


Measurements 


The  pressure  in  the  chamber  at  each  end  of  the 
rod  and  in  the  confining  pressure  chamber  sur¬ 
rounding  the  rod  were  measured  as  a  function  of 
time.  The  measurements  were  made  with  commer¬ 
cially  available  diaphragm-type  pressure  gages 
(Kulite  HKS-375). 

Examples  of  typical  transient  pressure  histories 
are  shown  in  Fig.  8.  Figures  8(a)  and  8(b)  show 
that  in  this  experiment  the  axial  unloading  oc¬ 
curred  in  about  25  fis  and  was  simultaneous  at  the 
two  ends.  The  radial  pressure,  shown  in  Fig.  8(c), 
was  constant  until  the  stress  waves  in  the  rod 
reached  the  gage  location  (midpoint  of  the  rod), 
and  remained  within  1  MPa  of  the  initial  value. 
The  variation  in  radial  pressure  appears  to  follow 
the  circumferential  strain  (not  shown),  suggesting 
that  the  pressure  variation  is  caused  by  volume 
changes  in  the  specimen,  produced  by  the  Poisson 
effect  as  the  axial  stress  wave  propagates  along  the 
rod.  This  result  shows  that  a  one-dimensional 
wave  analysis  is  not  rigorously  correct  for  these 
experiments  because  the  radial  stress  is  coupled  to 
the  axial  stress. 

Axial  and  circumferential  strains  were  mea¬ 
sured  at  several  axial  locations  on  the  surface  of 
the  rod,  using  commercially  available  2.5  cm-long 
foil-type  strain  gages  (Micromeasurements  MM- 
EA-06-10CBE-120).  This  length  of  strain  gage  re- 
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Fig.  8.  Pressure  measurements  in  experiment  3  (10  31  h\dro- 
static  preload). 


suits  in  considerable  averaging  of  strain  but.  for 
an  inhomogeneous  material  such  as  concrete,  the 
averaging  is  necessary.  Each  measurement  was 
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made  by  averaging  the  outputs  from  three  gages 
mounted  on  120  degree  intervals  at  the  same  axial 
position.  Averaging  these  three  measurements 
eliminates  bending  contributions  to  the  strain 
caused  by  small  curvature  in  the  rods.  Axial  strain 
was  measured  at  four  locations:  10  cm  from  each 
end  of  the  rod  and  7.6  cm  from  the  midpoint  on 
both  halves  of  the  rod.  Circumferential  strain  was 
measured  at  only  the  symmetric  points  7.6  cm 
from  the  midpoint. 

Concrete  rod  specimens 

The  concrete  tested  in  the  current  work  was 
made  from  graded  aggregate,  Portland  cement, 
and  water.  The  aggregate  was  local  river  rock  with 
rounded  shapes,  meeting  the  ASTM  C33  specifica¬ 
tion  for  size  distribution.  It  was  sieved  to  remove 
all  aggregates  not  passing  a  0.635  cm  (j  inch) 
opening.  The  average  static  compressive  strength 
was  about  60  MPa.  and  the  average  splitting  tensile 
strength  was  about  3.4  MPa.  The  average  elastic 
modulus  was  25  MPa,  and  Poisson's  ratio  was  0.2. 

The  rods  were  cast  about  90  cm  long,  and  then 
trimmed  to  76.2  cm.  Brass  sleeves,  measuring  2.5 
cm  long  and  0.04  cm  in  thickness,  were  epoxied  on 
to  the  rods  at  each  end.  These  sleeves  were  the 
sealing  surfaces  for  the  O-rings  in  the  testing 
apparatus. 

4.  Results 

The  four  experiments  described  in  this  paper 
are  listed  in  Table  1.  In  all  of  the  experiments  the 
strain  rate  in  the  front  of  the  rarefaction  waves 
was  about  10/s,  so  the  strain  rate  in  the  region  of 
superposition  of  the  waves  was  about  20/s.  Usu¬ 
ally  a  single  fracture  occurred  near  the  midpoint 
of  the  rod.  In  some  experiments  multiple  fractures 
occurred.  Based  on  measurements  of  axial  strains 
on  the  concrete  rods,  the  observed  strength  en¬ 
hancement  at  this  strain  rate  was  considerable. 

Experiments  on  unconfmed  rods 

Experiments  1  and  2  were  performed  on  con¬ 
crete  rods  with  no  confinement  and  a  static  axial 
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Fig.  9.  Axial  strain  records  at  -7.6  cm  (7.1  cm  from  fracture 
location)  in  experiment  1  (10.55  MPa  umaxjal  preload) 

preload  of  10.55  MPa.  Both  rods  failed  in  dy¬ 
namic  tension  at  a  single  location  within  0.5  cm 
from  the  midpoint,  and  no  secondary  damage  was 
visible. 

The  axial  strain  histories  from  Experiment  1 
measured  at  the  ±7.6  cm  locations  (referenced 
from  the  midpoint  of  the  rod)  are  shown  in  Figs.  9 
and  10.  Some  of  the  effects  of  iahomogeneity  in 
the  specimen  were  eliminated  by  scaling  the  re¬ 
corded  strain  signals  to  make  the  initial  strains 
(from  the  static  preload)  at  each  location  equal  to 
the  average  of  all  the  initial  strains.  On  the  time 
scale  of  these  plots,  the  explosive  charge  was  ini¬ 
tiated  at  t  =  0.100  ms.  Geometric  dispersion  in  the 
rarefaction  waves  as  they  propagate  from  the  ends 
of  the  rod  to  these  locations  causes  the  oscillation 
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Fig.  10.  Axtal  strain  records  at  7.6  cm  (8.1  cm  from  fracture 
location)  in  expenment  1  (10.55  MPa  umaxjal  preload). 


in  the  strain  records  just  before  the  arrival  of  the 
second  wave.  The  strain  rate  at  the  front  of  the 
rarefaction  waves  was  about  10/s,  so  the  strain 
rate  at  the  failure  location  was  about  20/s. 

The  axial  strains  show  the  effect  of  tensile 
failure  at  about  t  =  0.380  ms,  when  the  tensile 
strains  reach  a  peak  and  no  longer  follow  the 
history-  expected  for  elastic  response.  The  three 
individual  strains  at  both  locations  are  fairly  uni¬ 
form  even  after  the  effects  of  tensile  failure  arrive. 
After  fracture,  the  stress  waves  propagate  and 
reflect  in  the  two  separated  half-rods  with  stress- 
free  ends.  However,  reflections  do  not  return  to 
these  locations  during  the  time  period  shown  in 
these  figures. 

The  peak  average  strain  at  -7.6  cm  (7.1  cm 
from  the  failure  location)  was  160  microstrain. 
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Fig.  11.  Axial  strain  records  at  -7.6  cm  (7.5  cm  from  fracture 
location)  in  expenment  2  (10.55  MPa  uruaxiai  preload.) 


The  peak  average  strain  at  +7.6  cm  (8.1  cm  from 
the  failure  location)  was  210  microstrain.  Assum¬ 
ing  linear  elastic  response  at  the  gage  locations 
and  using  the  static  elastic  constants  for  the  con¬ 
crete,  the  higher  of  these  measurements  corre¬ 
sponds  to  an  axial  stress  of  about  5  MPa.  nearly 
50%  higher  than  the  average  static  splitting  tensile 
strength  for  this  family  of  rods. 

The  axial  strains  measured  in  Experiment  2  at 
±7.6  cm  are  shown  in  Figs.  11  and  12.  The 
individual  strains  show  wider  variations  than  those 
in  Experiment  1,  especially  after  the  effects  of 
tensile  failure  arrive.  However  the  average  strains 
are  very  nearly  the  same  as  those  in  the  previous 
experiment.  The  peak  average  strain  at  —7.6  cm 
(7.5  cm  from  the  failure  location)  was  170  micro- 
strain.  and  at  +7.6  cm  (7.7  cm  from  the  failure 
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would  be  high  in  both  cases  but  not  necessarily  by 
the  same  amount.) 

Experiments  on  confined  rods 

Experiments  3  and  4  were  performed  on  con¬ 
crete  rods  with  hydrostatic  preloads  of  10.31  MPa 
and  10.20  MPa,  respectively.  The  rod  in  Experi¬ 
ment  3  failed  in  dynamic  tension  at  two  locations. 
+  2.2  cm  and  — 15.2  cm.  No  large  voids  existed  at 
either  section.  The  rod  in  Experiment  4  failed  only 
at  -3.8  cm.  One  fairly  large  void  (0.8  cm  diam¬ 
eter)  was  evident  at  the  failed  section. 

Plots  of  the  average  axial  strains  at  ±  7.6  cm  in 
Experiment  3  are  shown  in  Fig.  13.  The  axial 
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Fig  12.  Axial  strain  records  at  7.6  cm  (7.7  cm  from  fracture 
location)  in  experiment  2  (10.55  MPa  uniaxial  preload.) 


location)  it  was  190  microstrain.  The  elastic  axial 
stress  corresponding  to  the  highest  measured  strain 
is  about  4.7  MPa,  about  40%  higher  than  the 
average  static  splitting  tensile  strength  for  this 
family  of  rods. 

For  a  material  like  concrete,  the  waveforms  and 
peak  strains  measured  in  Experiments  1  and  2 
demonstrate  very  good  reproduciblity  of  the  re¬ 
sults.  It  is  recognized,  however,  that  the  40-50% 
estimated  strength  enhancement  is  much  less  than 
the  400%  enhancement  reported  by  Birkimer 
(1968)  for  these  strain  rates.  This  contrast  may 
simply  reflect  the  differences  between  the  materi¬ 
als  tested.  It  may  also  be  associated  with  the 
assumption  of  elastic  strains.  (If  the  measured 
strains  were  inelastic,  the  estimates  of  strength 


Fig.  13  Average  axial  strains  at  +7  6  cm  in  experiment  3 
(10.31  MPa  hydrostatic  preload) 
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Table  1 

Dynamic  tension  experiments  on  concrete  rods 


i:3 


Experiment 
number  1 

AjuaJ  preload 
(MPa) 

Radial  preload 
(MPa) 

Elastic  modulus  2 
(GPa) 

Dynamic  tensile 
strength  *  (MPa) 

Fracture  locations 
(cm  from  midpoint) 

i 

10.55 

0.0 

24.1 

5.0 

+  0.46 

2 

10.55 

0.0 

24.7 

4.7 

+  0.10 

3 

10.31 

10.31 

2fT.4 

6.9 

-15.2.  +2.2 

4 

10.20 

10.20 

23.0 

6.7 

-3.8 

1  These  experiments  were  originally  titled  Test  42.  43,  44.  and  46,  respectively  (Gran.  1985). 

2  The  elastic  constants  were  determined  from  the  static  preload  and  the  average  initial  strains. 

3  These  are  estimates  corresponding  to  the  elastic  stresses  computed  from  the  peak  measured  tensile  strains. 


strain  corresponding  to  zero  axial  stress  is  145 
microstrain.  The  failure  at  +  2.2  cm  produced  the 
drop  in  axial  strains  at  ±  7.6  cm  at  about  t  =  0.380 
ms.  The  fracture  at  - 15.2  cm  produced  a  new  free 
surface  there,  and  the  rarefaction  wave  reflected 
from  this  section  and  returned  to  the  strain  gages 
at  -7.6  cm.  The  other  half  of  the  rod  remained 
intact,  so  no  reflections  from  the  end  of  the  rod 
returned  during  this  period. 

The  peak  average  strain  at  —7.6  cm  was  390 
microstrain.  At  +7.6  cm  it  was  315  microstrain. 
Again  assuming  linear  elastic  response  at  the  gage 
locations,  using  the  elastic  modulus  given  in  Table 
1,  and  accounting  for  the  nominal  radial  stress  of 
10.20  MPa,  these  axial  strains  correspond  to  axial 
stresses  of  about  6.9  MPa  and  5.6  MPa,  respec¬ 
tively.  The  higher  of  these  stresses  is  about  100$ 
greater  than  the  average  static  splitting  tensile 
strength  for  this  family  of  rods.  It  is  also  about 
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Fig  14.  Average  axial  strain  at  -  7.6  cm  (3.8  cm  from  fracture) 
in  experiment  4  (10.31  MPa  hydrostatic  preload.) 


40$  higher  than  the  unconfined  tensile  strength  at 
the  same  strain  rate,  observed  in  Experiments  1 
and  2. 

In  Experiment  4,  the  unloading  device  at  one 
end  did  not  perform  properly.  The  pressure 
dropped  smoothly  in  about  30  (is,  but  it  remained 
at  zero  for  only  about  50  (is  before  recompression 
occurred.  Posttest  inspection  of  the  unloading  de¬ 
vice  at  this  end  showed  that  the  segmented  ring 
was  only  partially  removed  by  the  explosive  and 
was  trapped  between  the  support  tube  and  the 
spacer.  Apparently,  this  limited  the  travel  of  the 
piston  so  that  the  extension  of  the  rod  recom¬ 
pressed  the  oil.  In  addition  to  this  problem,  the 
strain  gages  at  +  7.6  cm  did  not  function. 

The  axial  strain  measured  at  -7.6  cm  in  Ex¬ 
periment  4  is  shown  in  Fig.  14.  The  wavefront  of 
the  first  rarefaction  shows  a  longer  rise  time  than 
was  typical  of  the  waves  in  previous  experiments, 
possibly  because  of  the  malfunction  of  the  unload¬ 
ing  device.  However,  after  the  arrival  of  the  sec¬ 
ond  rarefaction  wave  (beginning  at  about  r  =  0.365 
ms),  the  tensile  strain  rate  is  about  the  same  as  in 
Experiment  3.  The  peak  axial  strain  is  470  micro¬ 
strain,  corresponding  to  an  elastic  axial  stress  of 
6.7  MPa.  This  is  about  the  same  as  that  computed 
from  the  strains  in  Experiment  3  (the  elastic  mod¬ 
uli  were  significantly  different). 

Thus,  although  there  were  difficulties  with 
multiple  fracture,  malfunctions  of  the  unloading 
device,  and  loss  of  some  strain  records.  Experi¬ 
ments  3  and  4  demonstrated  the  capability  of  the 
experimental  technique  to  produce  tensile  failure 
at  a  strain  rate  of  10  to  20/s  with  independently 
controlled  confining  pressure.  The  surprising  re¬ 
sult  that  the  apparent  dynamic  tensile  strength  is 
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enhanced  by  the  confinement  is  in  contrast  to 
static  data,  e.g.  (Saucier,  1974),  but  no  other  dy¬ 
namic  data  are  available  for  comparison. 

5.  Interpretation  of  the  results 

Assuming  the  peak  observed  strains  in  these 
experiments  to  be  elastic,  the  unconfined  tensile 
strength  of  the  concrete  at  a  strain  rate  of  10  to 
20/s  averaged  over  40%  higher  than  the  static 
splitting  tensile  strength.  At  the  same  strain  rate, 
the  tensile  strength  with  10  MPa  confining  pres¬ 
sure  was,  on  average,  about  100%  higher  than  the 
static  splitting  tensile  strength,  and  about  40% 
higher  than  the  unconfined  tensile  strength  at  10 
to  20/s. 

These  strength  estimates  are  based  on  the  as¬ 
sumption  that  the  measured  strains  are  elastic,  but 
this  assumption  may  not  be  justified.  Thus,  to 
further  interpret  the  experimental  results  wave- 
propagation  calculations  were  performed  using  a 
one-dimensional  strain-softening  model  (Gran, 
1985).  The  model  is  based  on  the  assumption  that 
the  stress-strain  relation  is  not  a  property  of  a 
material  point,  but  is  an  average  property  of  a 
finite  volume  of  material  containing  a  developing 
crack  or  failure  plane.  The  stress-strain  relation 
thus  has  associated  with  it  a  finite  dimension, 
namely  the  average  crack  separation  distance. 
Using  this  model,  the  two  dynamic  unconfined 
tension  experiments  were  simulated,  and  by  trial- 
and-error  good  agreement  with  the  measured 
strains  was  obtained  in  both  cases. 

Whereas  the  dynamic  experiments  produced  a 
single  fracture  plane  with  no  visible  secondary 
damage,  the  calculations  predicted  some  inelastic 
strain  to  occur  virtually  thoughout  the  specimen, 
with  concentrations  of  inelastic  strain  within  a  few 
centimeters  of  the  locations  of  complete  fracture. 
In  addition,  the  calculations  suggest  that  the  strain 
history  measured  a  few  centimeters  from  the  loca¬ 
tion  of  fracture  is  primarily  a  function  of  inelastic 
wave  propagation  from  the  fracture  location  to  the 
strain  gage  (through  the  sites  of  concentrated  in¬ 
elastic  strain),  and  is  less  dependent  on  the  behav¬ 
ior  of  the  material  right  at  the  fracture.  However, 
they  also  showed  that  the  strains  were  measured  at 


locations  where  the  inelasticity  was  slight,  so  that 
in  this  case  the  estimates  of  strength  are  reasona¬ 
ble.  The  nonlinear  analyses  performed  by  Gran 
(1985)  will  be  described  in  a  subsequent  paper. 

Even  without  nonlinear  analyses,  however,  these 
experiments  show  that  the  assumption  of  elastic 
response  everywhere  except  at  the  failure  location 
is  not  valid,  even  when  only  one  section  was 
visibly  damaged.  According  to  that  assumption,  as 
depicted  in  Fig.  4,  the  observed  peak  strain  should 
be  £(i  +  a)  times  the  prestrain,  where  a  is  the 
ratio  of  the  strength  to  the  preload.  That  is.  the 
observed  peak  strain  should  never  be  less  than 
half  the  prestrain,  even  when  the  strength-to-pre- 
load  ratio  is  zero.  (The  measurements  need  to  be 
made  at  distances  greater  than  JA(1  -  a)  from  the 
failure  location  in  order  to  not  be  limited  by 
residual  pre-compression.)  The  observed  peak 
strains  in  the  unconfined  experiment  were  all  less 
than  half  the  prestrain,  and  the  measurement  loca¬ 
tions  were  well  outside  the  region  where  the  strain 
would  be  limited  by  residual  precompression. 

Similarly,  in  the  confined  experiments  the  peak 
observed  strains  at  ±  7.6  cm  would  correspond  to 
stresses  equal  to  ](1  +  a)a0.  For  peak  stresses  of. 
say,  6.8  MPa  (computed  from  the  peak  measured 
strains)  and  an  intitial  preload  of  10  MPa.  a 
would  be  0.36.  The  strength  at  the  failure  loca¬ 
tion,  given  by  ao0,  would  only  be  3.6  MPa.  or 
about  half  the  peak  stress  occurring  remote  from 
the  failure  location  as  estimated  from  measured 
strains. 

Thus,  estimates  of  tensile  strength  based  on 
elastic  wave  analyses  are.  in  general,  not  valid  for 
this  type  of  experiment. 

6.  Summary  and  conclusions 

An  experimental  technique  to  measure  the 
tensile  response  of  brittle  geologic  sobds  at  strain 
rates  of  approximately  10  to  20/s  was  developed 
and  applied  to  concrete  rod  specimens.  In  all  of 
the  experiments  the  primary  failure  location  oc¬ 
curred  within  a  few  centimeters  of  the  midpoint  of 
the  rod.  In  some  of  the  experiments  a  second 
failure  occurred  near  one  quarter-point.  Successful 
measurements  were  made  of  the  applied  pressures 
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and  the  surface  strains.  The  effects  of  tensile  failure 
on  the  stress  waves  in  the  specimen  were  observed 
in  the  surface  strain  measurements. 

The  initial  results  appear  to  be  sufficiently  re¬ 
producible  to  warrant  a  fairly  detailed  interpreta¬ 
tion.  The  input  waveforms,  measured  at  ±  7.6  cm 
from  the  midpoint  of  the  rod,  are  nearly  identical 
from  test  to  test.  For  similar  initial  conditions,  the 
average  axial  strain  measurements  after  fracture 
are  also  very  similar  in  form  and  magnitude. 

In  general,  the  experimental  results  require  a 
nonlinear  analysis  for  proper  interpretation.  How¬ 
ever,  for  the  particular  experiments  described  here, 
a  simple  elastic  interpretation  of  peak  strains  gives 
reasonable  estimates  for  tensile  strength  because 
the  inelastic  strains  at  the  measurement  locations 
are  not  large.  This  simple  interpretation  indicates 
that  there  is  about  a  40 %  increase  in  the  tensile 
strength  of  this  particular  concrete  at  a  strain  rate 
of  10  to  20/s,  compared  to  the  static  splitting 
strength.  The  data  ..Iso  indicate  that  there  is  an 
additional  40%  increase  of  the  dynamic  tensile 
strength  with  a  static  confining  pressure  of  10 
MPa.  However,  these  estimates  are  only  ap¬ 
proximate  because  inelastic  response  in  the  speci¬ 
men  is  not  included  in  this  interpretation. 

Clearly,  more  experiments  and  interpretive 
analyses  are  needed  to  even  empirically  char¬ 
acterize  the  tensile  failure  of  concrete  at  these 
strain  rates.  In  addition,  the  application  of  the 
experimental  method  to  brittle  geologic  materials, 
particularly  fine-grained  rocks,  is  an  obvious  fu¬ 
ture  direction.  Coupled  with  nonlinear  wave  prop¬ 
agation  analyses,  this  experimental  method  should 
provide  unique  and  valuable  data  for  the  develop¬ 
ment  of  rate-dependent  tensile  failure  models  for 
these  materials. 
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STRAIN-SOFTENING  CALCULATIONS  FOR  CONCRETE 
IN  DYNAMIC  UNIAXIAL  TENSION 


By  James  K.  Gran1  and  Lynn  Seaman2,  M.ASCE 


ABSTRACT:  A  one-dimensional  strain-softening  model  was  used  in  wave-propagation 
calculations  to  interpret  dynamic  unconfined  tension  experiments  on  concrete  rods.  The 
model  is  based  on  the  assumption  that  the  stress-strain  relation  is  not  a  property  of  a 
material  point,  but  is  an  average  property  of  a  finite  volume  of  material  containing  a 
developing  crack  or  failure  plane.  The  stress-strain  relation  has  associated  with  it  a  finite 
dimension,  namely  an  effective  crack  separation  distance.  Two  experiments  were  simu¬ 
lated,  and  with  suitable  choice  of  the  material  parameters  good  agreement  with  the  meas¬ 
ured  axial  strain  histories  was  obtained  in  both  cases.  In  addition  to  providing  an  estimate 
of  the  dynamic  tensile  properties  of  the  concrete,  these  calculations  suggest  that  tensile 
damage  in  the  concrete  was  distributed  over  several  centimeters  surrounding  the  location  of 
fracture.  The  strain  history  measured  a  few  centimeters  from  the  fracture  appears  to  be 
primarily  a  function  of  inelastic  wave  propagation  from  the  fracture  to  the  strain  gage,  and 
is  only  weakly  dependent  on  the  behavior  of  the  material  at  the  location  of  complete  frac¬ 
ture. 

SUMMARY:  A  one-dimensional  strain-softening  model  was  used  in  wave-propagation 
calculations  to  interpret  results  of  dynamic  tension  experiments  on  concrete  rods.  Two 
dynamic  unconfined  tension  experiments  were  simulated,  and  good  agreement  with  the 
measured  axial  strain  histories  was  obtained  in  both  cases. 

KEY  WORDS:  strain-softening,  calculations,  concrete,  dynamic,  tension,  fracture,  crack¬ 
ing,  damage 
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INTRODUCTION 


Tensile  failure  in  concrete  and  brittle  geologic  materials  results  from  the  growth  and 
coalesence  of  tensile  microcracks  with  a  corresponding  loss  of  tensile  strength.  The  under¬ 
standing  and  quantification  of  this  failure  process  are  currently  not  sufficient  to  provide  a 
detailed  micromechanical  model.  This  situation  prompts  the  use  of  so-called  continuum 
strain-softening  models  in  numerical  simulations  of  tensile  failure.  These  models  allow 
decreasing  stress  with  increasing  strain  by  averaging  over  a  finite  region  the  deformation 
associated  with  cracking. 

Some  notable  examples  of  continuum  strain-softening  models  that  have  been  utilized 
successfully  to  describe  quasi-static  experimental  behavior  are  the  line  crack  model  intro¬ 
duced  by  Hillerborg,  et  al  (1976),  the  crack  band  model  of  Bazant  and  Oh  (1983),  and  the 
composite  fracture  model  developed  by  Wiliam,  et  al  (1984).  As  pointed  out  by  Ottosen 
(1986),  each  of  these  models  has  the  necessary  feature  that  a  fundamental  relationship 
between  stress  and  displacement  across  a  softening  zone  is  preserved  as  the  size  of  the 
numerical  cell  containing  the  softening  zone  is  varied.  Bazant  (1976)  concluded  that  the 
proper  size  of  the  softening  cell  must  be  determined  experimentally. 

Recently,  the  proper  construction  and  use  of  strain-softening  models  for  dynamic 
problems  has  been  a  topic  of  considerable  interest.  Sandler  and  Wright  (1984)  demon¬ 
strated  with  numerical  and  analytical  examples  that  any  attempt  to  model  uniform  deforma¬ 
tion  in  a  rate-independent  strain-softening  continuum  will  always  result  in  an  unbounded 
strain  concentration.  Belytschko  and  Bazant  (1984)  also  used  numerical  and  analytical 
examples  to  show  that  for  a  given  stress-strain  relation  with  softening,  the  calculated  work 
done  to  fail  the  material  is  cell-size  dependent.  Furthermore,  as  the  cell  size  tends  to  zero, 
the  work  also  tends  to  zero,  which  is  certainly  not  a  physically  plausible  result.  Both  of 
these  studies  considered  only  softening  continuua,  that  is,  rather  than  preserving  a  stress- 
displacement  relationship  across  a  softening  zone,  these  models  preserve  a  stress-strain 
relationship  as  the  size  of  the  cell  is  varied.  A  conclusion  of  both  studies  is  that  such  a 
model  is  physically  implausible. 

This  paper  describes  an  attempt  to  apply  a  one-dimensional  model  of  the  Hillerborg- 
Bazant-Willam  type  to  a  dynamic  problem,  and  to  assess  the  ability  of  the  model  to  repro¬ 
duce  dynamic  experimental  observations.  We  assume  that  the  physical  source  of  strain 
softening  is  tensile  cracking,  that  the  deformation  is  naturally  concentrated  at  the  cracks, 
and  that  the  stress  versus  crack-opening  relation  is  a  fundamental  property  of  the  material. 
In  addition,  we  assume  that  the  softening  stress-strain  relation  is  not  a  property  of  a 
material  point,  but  is  an  average  property  of  a  finite  volume  of  material  containing  a 
developing  crack  or  failure  plane.  The  stress-strain  relation  thus  has  associated  with  it  a 
characteristic  finite  dimension,  namely  the  effective  crack  separation  distance. 


The  strain-softening  calculations  described  here  were  guided  by  experiments  per¬ 
formed  by  Gran  et  al  (1987),  whose  experimental  technique  is  illustrated  in  Figure  1.  With 


this  technique  a  cylindrical  concrete  rod  specimen  is  first  loaded  in  static  uniaxial  compres¬ 
sion,  then  the  axial  load  is  released  from  each  end  simultaneously  and  very  rapidly.  The 
resulting  relief  waves  interact  in  the  center  of  the  rod  to  produce  a  dynamic  tensile  stress 
equal  in  magnitude  to  the  original  static  compression  at  a  strain  rate  of  10/s  to  20/s.  Typi¬ 
cally  in  these  experiments,  a  single  complete  fracture  across  the  rod  occurs  near  the  mid¬ 
point  of  the  rod.  Transient  measurements  are  made  of  the  axial  load  at  each  end  of  the  rod 
and  of  axial  and  circumferential  surface  strains  at  several  locations  along  the  length  of  the 
rod.  An  example  of  axial  strain  measurements  from  a  typical  experiment  is  shown  in  Fig¬ 
ure  2.  The  static  compressive  preload  in  this  experiment  was  10  MPa;  the  nominal  static 
tensile  strength  of  the  concrete  was  3.4  MPa.  The  strains  shown  in  Figure  2  were  meas¬ 
ured  at  symmetric  locations  7.62  cm  from  the  midpoint  of  the  rod  using  three  2.54-cm-long 
strain  gages  at  each  location.  A  single  complete  fracture  across  the  rod  occurred  at  0.1  cm 
from  the  midpoint  of  the  rod,  but  no  other  macroscopic  tensile  damage  was  visible. 

These  experimental  results  do  not  give  directly  the  strength  of  the  concrete  or  the 
strength  reduction  as  a  function  of  accumulating  tensile  damage.  Thus,  the  objective  of  the 
calculations  described  in  this  paper  is  to  estimate  for  these  experiments  the  tensile  strength 
and  the  strength  reduction  as  a  function  of  accumulating  damage  exhibited  by  pertinent 
regions  of  the  concrete  specimens.  The  analyses  were  performed  by  numerically  integrat¬ 
ing  the  governing  equations  for  stress  waves  in  the  rod,  using  an  assumed  strain-softening 
response.  By  comparing  the  computed  strain  histories  with  the  measurements  from  the 
experiments,  the  assumed  strain-softening  response  was  adjusted  to  give  good  correlation 
with  the  observations.  Although  this  exercise  does  not  produce  a  unique  constitutive 
model  for  concrete  tensile  failure,  it  does  provide  insight  into  the  dynamic  tensile  failure  of 
concrete  by  making  an  estimate  of  what  must  have  happened  in  the  neighborhood  of  the 
complete  fracture  location-a  part  of  the  response  that  cannot  be  measured  directly. 

FORMULATION  OF  THE  STRAIN-SOFTENING  MODEL 

Our  approach  is  to  assume  that  the  rod  comprises  a  single  row  of  one-dimensional 
material  cells,  each  containing  a  single  potential  fracture  plane  (crack).  The  analytical 
model  used  for  tensile  failure  is  simply  a  prescribed  relation  between  the  stress  in  a  cell 
and  the  fracture  volume  per  unit  area  (crack  opening)  created  during  the  failure  process.  A 
pictorial  representation  of  the  assumed  failure  process  is  shown  in  Figure  3.  The  material 
is  initially  elastic  with  modulus  E,  that  of  intact  material.  Fracture  volume  per  unit  area 
(5)  is  created  only  after  the  stress  reaches  the  tensile  strength  (Co).  Then  under  continued 
tensile  strain,  the  fracture  volume  per  unit  area  increases  to  a  critical  value  (8C)  as  the 
strength  is  reduced  to  zero,  resulting  in  complete  fracture.  Note  that  5  has  the  units  of  dis¬ 
placement  and  can  be  regarded  as  an  effective  crack  opening.  The  relation  between  C  and 
5  (for  loading  and  unloading)  is  a  basic  property. 

The  initial  length  of  the  material  cell  is  Iq.  This  is  a  characteristic  dimension  of  the 
material  and  can  be  regarded  as  an  effective  crack  spacing  in  this  one-dimensional 
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(c)  Tensile  failure  at  6  =  Sc. 
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Figure  3.  Concept  of  the  one-dimensional  tensile  failure  model 


B-7 


formulation.  The  source  of  the  effective  crack  spacing  may  be  the  distribution  of  aggregate 
grains  as  crack  nucleation  sites.  Another  source  of  effective  crack  spacing  could  be 
stress-wave  interactions  that  permit  the  growth  of  only  those  cracks  that  are  sufficiently 
separated,  as  suggested  by  the  work  of  Kipp  and  Grady  (1985). 

An  example  of  a  relation  between  a  and  8  is  shown  in  Figure  4(a)  as  a(5).  In  this 
example,  the  strength  shows  no  increase  beyond  the  initial  strength  and  decays  to  zero 
linearly.  Unloading  and  reloading  follow  a  straight  line  leading  to  and  from  the  origin, 
implying  that  some  of  the  work  expended  in  creating  the  fracture  is  recoverable.  As 
shown  in  Figure  4(b),  the  strain  value  corresponding  to  complete  separation  depends  on  the 
original  dimensions  of  the  fracturing  cell.  That  is,  for  a  given  relation  between  stress  and 
fracture  volume  per  unit  area,  the  relation  between  stress  and  strain  depends  on  the  crack 
spacing  (cell  size)  assumed  in  the  numerical  discretization. 

As  mentioned  above,  when  implemented  as  a  stress-strain  relation,  this  model  is 
essentially  the  same  as  those  described  by  Hillerborg,  et  al  (1976),  Bazant  (1976),  and  Wil¬ 
iam,  et  al  (1984).  In  the  softening  regime  (e>0(/E),  the  fracture  volume  per  unit  area  is 
8  =  (e-a/E)/0.  In  general,  the  stress  reduction  is  given  by  0(8),  and  the  stress-strain  rela¬ 
tion  is  obtained  by  inserting  the  expression  for  8  and  solving  for  O(e).  For  the  special 
case  of  linear  stress  reduction  shown  in  Figure  4  the  softening  stress-strain  relation  is 

£o 

o=E- - —  (Cuu-e)  (i) 

which  depends  on  /0  through  the  definition  of  =  Sc//q-  An  explicit  expression  for  8  in 
terms  of  e  is  then 

8=/0 — — —  (e-Eo)  (2) 

EulfEo 

To  summarize,  the  tensile  failure  model  is  defined  by  the  elastic  modulus  E,  the  strength 
Oq,  the  relation  between  stress  and  fracture  volume  per  unit  area  0(8),  and  the  effective 
crack  spacing  /0. 

Even  though  the  stress-strain  relation  depends  on  the  cell  size,  the  work  done  per  unit 
area  to  completely  fracture  a  cell  is  independent  of  the  cell  size  because  all  the  work  goes 
into  creating  the  fracture,  and  the  elastic  material  returns  to  a  state  of  zero  internal  energy. 
For  linear  softening  the  work  done  per  unit  area  to  completely  fracture  a  cell  is 
Gc  =  O08c/2.  Thus,  the  work  done  per  unit  area  to  completely  fracture  a  cell  depends  only 
on  the  critical  stress  and  the  critical  fracture  volume  per  unit  area,  and  does  not  depend  on 
the  cell  size. 


However,  if  a  cell  is  only  partially  damaged,  that  is  strained  to  Ej^,  where  Eniax<eult» 
then  the  unrecovered  work  per  unit  area  for  linear  unloading  to  the  origin  is 


G'  = 


1  _  c  (Emax- ^0) 
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Figure  4.  Implementation  of  the  tensile  failure  model  as  a  strain-softening 


which  depends  on  Iq  through  e^.  This  dependence  occurs  because  for  a  given  maximum 
strain  the  maximum  fracture  volume  per  unit  area  (crack  opening)  depends  on  the  cell  size 
(crack  spacing),  as  given  in  equation  (2). 

It  should  be  expected,  therefore,  that  strains  calculated  with  this  model  in  wave- 
propagation  analyses  will  depend  on  the  cell  size  (crack  spacing)  chosen  for  the  calcula¬ 
tion.  Indeed,  because  the  cell  size  represents  an  allowed  crack  spacing,  it  is  a  material 
parameter  that  must  be  determined.  This  crack  spacing  is  similar  to  the  crack  band  width 
wc,  considered  by  Bazant  and  Oh  (1983).  Consequently,  as  part  of  the  parameter  study 
with  the  tensile  failure  model,  a  study  was  made  of  cell  size  dependence.  Ottosen’s  (1986) 
restriction  on  cell  size,  namely  that  Iq< 2GcE/Cq  .  is  satisfied  by  requiring  that  Eui^Eo- 

Finally,  we  point  out  that  the  model  for  the  failure  process  used  here  is  not  rate 
dependent.  Although  we  believe  the  true  response  to  be  rate  dependent,  the  range  of  strain 
rates  in  the  experiments  under  study  is  not  more  than  a  factor  of  2  (10/s  to  20/s).  Thus,  as 
a  simplification,  rate  dependence  has  not  been  included  explicitly,  although  the  parameters 
chosen  by  trial  and  error  to  match  the  experimental  results  define  our  estimate  of  the 
material  behavior  at  strain  rates  of  10/s  to  20/s. 


PARAMETER  STUDY 

Several  one-dimensional  stress  wave-propagation  calculations  were  performed  using  a 
finite-difference  scheme  (Seaman,  1978)  to  integrate  the  governing  equations.  These  calcu¬ 
lations  illustrate  the  sensitivity  of  the  numerical  results  to  the  parameters  defining  the 
assumed  strain-softening  properties.  In  one  set  of  calculations,  the  properties  of  the 
material  at  the  location  of  complete  fracture  were  varied  while  the  properties  in  the  rest  of 
the  rod  remained  fixed.  This  set  of  calculations  establishes  requirements  for  controlling  the 
location  of  complete  fracture.  In  another  set  of  calculations,  the  unloading  modulus  was 
varied  to  determine  its  effect  on  calculated  strain  histories.  In  a  third  set  of  calculations, 
the  cell  size  was  varied  as  the  0(5)  relation  remained  fixed.  These  sets  of  calculations 
demonstrate  the  qualitative  nature  of  the  model. 

The  correlation  of  the  predictions  with  the  experimental  results  was  made  by  compar¬ 
ing  the  calculated  axial  strains  with  the  measurments  made  in  an  experiment  by  Gran,  et  al 
(1987).  The  strain  histories  from  the  selected  experiment  were  shown  in  Figure  2.  The 
peak  strain  (measured  7.7  cm  from  the  location  of  complete  fracture)  was  190  microstrain, 
which,  if  elastic,  would  correspond  to  an  axial  stress  of  about  4.7  MPa  (40%  greater  than 
the  static  strength).  In  each  calculation,  a  static  preload  of  10.20  MPa  was  applied  at  each 
end  of  the  rod,  and  removed  with  an  exponential  decay  with  a  time  constant  of  15  (is. 

Variation  of  the  Properties  at  the  Location  of  Complete  Fracture 

In  Calculation  A  every  cell  in  the  rod  was  given  the  property  shown  in  Figure  5(a) 


labeled  "Softening".  The  cell  width  was  0.333  cm  and  the  fracture  volume  per  unit  area  at 
complete  separation  was  about  7fim.  The  predicted  stress  and  strain  histories  for  this  cal¬ 
culation  are  not  shown  because  a  single  complete  fracture  at  the  origin  did  not  occur,  mak¬ 
ing  a  comparison  with  the  experimental  results  meaningless.  Complete  fracture  was  com¬ 
puted  to  occur  at  four  locations,  0.66  cm  and  3.33  cm  from  the  origin  in  both  halves  of  the 
rod.  A  segment  of  the  rod  about  4  cm  long  had  a  fairly  uniform  stress  distribution  at  the 
time  the  assumed  strength  was  first  reached.  Apparently,  numerical  perturbations  in  the 
unloading  waveform  caused  the  stress  to  first  exceed  the  strength  at  a  point  other  than  the 
origin. 

This  result  demonstrates  that  if  all  the  cells  have  the  same  strength  and  the  same 
strain-softening  relation,  the  calculated  location  of  complete  fracture  is  determined  by 
numerical  perturbations  in  the  load.  It  also  implies  that  experiments  in  which  the  loading  is 
presumably  symmetric,  but  in  which  complete  fracture  occurs  at  an  asymmetric  location  (a 
likely  occurrence)  cannot  be  correctly  represented  with  uniform  properties  along  the  length 
of  the  rod.  However,  we  know  that  the  properties  of  the  concrete  are  non-uniform.  Con¬ 
sequently,  in  the  experiment,  the  location  of  complete  fracture  is  controlled  by  the 
existence  of  a  locally  weak  section  near  the  midpoint.  Unloading  from  this  fracture 
quenches  the  growth  of  cracks  at  other  potential  fracture  locations.  (This  would  also  be 
true  in  a  static  experiment,  that  is,  the  weakest  section  in  the  specimen  would  control  the 
location  of  failure.  All  other  locations  would  unload  at  the  time  of  failure.)  Thus,  we 
adopt  the  standard  practice  of  introducing  a  perturbation  to  the  specimen  to  induce  strain 
localization  at  the  desired  location. 

Calculation  B  was  performed  to  demonstrate  how  the  properties  at  the  origin  can  be 
modified  to  ensure  that  complete  fracture  occurs  only  at  the  origin.  The  rod  was  modeled 
with  the  "Softening"  relation  used  before,  but  the  cell  at  the  origin  was  given  80%  of  the 
strength  of  the  cells  in  the  rest  of  the  rod,  using  the  strain-softening  relation  labeled  "80% 
Strength"  in  Figure  5(b).  In  this  calculation  complete  fracture  occurred  only  at  the  origin. 
The  strain  history  from  this  calculation  is  plotted  in  Figure  6.  The  shape  of  the  calculated 
strain  history  at  the  strain  gage  location  is  noticeably  flatter  than  the  measured  strain  his¬ 
tories.  Although  complete  fracture  occurred  only  at  the  origin,  the  plot  of  peak  strains  in 
Figure  7  shows  that  inelastic  tensile  strain  (associated  with  the  peak  computed  fracture 
volume  per  unit  area)  was  distributed  over  several  centimeters  and  exhibited  a  concentra¬ 
tion  at  the  2.66-cm  location. 

Variation  of  the  Unloading  Modulus 

Calculation  C  was  performed  to  investigate  the  effects  of  the  unloading  modulus 
assumed  for  the  cells  that  undergo  inelastic  strain  but  do  not  fracture  completely.  In  the 
previous  two  calculations,  unloading  followed  a  straight  line  to  the  origin.  In  Calculation 
C,  unloading  followed  a  straight  line  parallel  to  the  elastic  loading  line.  The  relations  used 
in  Calculation  C  are  plotted  in  Figure  5(c).  The  strain  histories  calculated  with  elastic 
unloading  are  shown  in  Figure  8.  These  should  be  compared  with  the  strain  histories  from 
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Calculation  B  (Figure  6).  Primarily,  the  effect  of  elastic  unloading  is  to  produce  a  residual 
tensile  strain.  This  puts  the  calculated  strain  at  7.62  cm  in  slightly  better  agreement  with 
one  of  the  measurements  from  the  experiment.  Another  effect  of  elastic  unloading  is  to 
decrease  the  duration  of  the  strain  pulse.  This  effect  is  caused  by  the  greater  wave  speed 
associated  with  the  steeper  unloading  slope. 

Variation  of  the  Cell  Size 

The  cell  size  in  these  calculations  has  the  role  of  a  material  parameter  because  it  is 
equivalent  to  an  effective  crack  spacing.  That  is,  each  cell  represents  a  potential  failure 
plane  surrounded  by  elastic  material.  In  the  derivation  of  the  strain-softening  model,  it  was 
pointed  out  that  the  ultimate  strain  consistent  with  the  critical  fracture  volume  per  unit  area 
depends  on  the  cell  size  in  the  calculation  £uit= 8c//q.  In  addition,  the  unrecovered  work 
per  unit  area  in  a  partially  damaged  cell  depends  on  the  cell  size  (Eq.  3). 

Thus,  before  the  simulations  of  the  experiments  were  attempted,  a  study  of  cell  size 
dependence  was  conducted.  Calculation  set  D  was  performed  with  cell  sizes  of  0.3  cm, 
0.6  cm,  and  1.0  cm.  The  same  0(8)  relation  used  in  all  the  calculations  was  linear  with 
Oc  =  5  MPa  and  Sc  =  6  (Im.  As  suggested  by  the  previous  results  of  the  parameter  study, 
the  cell  at  the  origin  was  given  80%  of  the  strength  in  the  rest  of  the  rod.  As  before,  the 
rod  was  preloaded  with  10  MPa  compression,  and  the  load  was  released  with  an  exponen¬ 
tial  decay  with  a  time  constant  of  15  {is.  The  calculated  strains  at  2  cm  from  the  origin 
are  plotted  in  Figure  9.  The  strains  from  the  0.3-cm,  and  0.6-cm  cells  are  in  good  agree¬ 
ment  at  both  locations  and  appear  to  approximate  the  ideal  plasticity  solution.  The  strains 
from  the  1-cm  cells  are  qualitatively  different  at  both  locations.  This  calculation  also 
predicts  a  very  large  local  strain  at  the  4-cm  location  that  is  not  predicted  in  the  other 
cases.  The  contrast  indicates  how  the  predicted  response  depends  on  the  effective  crack 
spacing  (cell  size). 

SIMULATIONS  OF  THE  EXPERIMENTS 

The  experiments  analyzed  here  are  labeled  Experiment  1  and  Experiment  2.  To  simu¬ 
late  an  experiment,  the  stress-strain  relations  and  the  cell  size  in  the  rod  and  at  the  loca¬ 
tions  of  complete  fracture  were  varied  until  the  calculated  strains  matched  the  measured 
strains.  Both  of  the  experiments  were  simulated  very  well  using  the  same  failure  proper¬ 
ties.  The  measured  elastic  constants  and  the  measured  density  were  used  in  both  cases. 
Complete  fracture  was  forced  to  occur  at  the  observed  locations  of  failure  by  degrading  the 
properties  of  the  rod  at  those  locations  to  80%  of  the  strength  in  the  rest  of  the  rod,  as 
described  below. 

The  simulations  shown  below  were  performed  using  0.635-cm  cells.  This  cell  size 
was  found  to  be  the  one  that  worked  best.  When  0.333-cm  cells  were  used,  the  roundness 
of  the  peaks  of  the  measured  strain  histories  could  not  be  matched.  The  calculated  strain 
histories  were  too  flat.  When  1-cm  cells  were  used,  additional  complete  fractures  would 


occur  in  the  rod.  The  0.635-cm  dimension  is  the  diamter  of  the  largest  aggregate  particles 
and  is  roughly  equal  to  the  average  axial  spacing  of  the  largest  aggregate  particles.  This 
cell  size  is  considerably  smaller  than  was  found  to  work  best  in  static  analyses  performed 
by  Bazant  and  Oh  (1983).  Perhaps  this  is  because  stress-wave  loading  excites  more  poten¬ 
tial  fracture  sites  before  the  unloading  wave  from  a  complete  fracture  quenches  further 
damage. 

The  strain-softening  relations  used  in  the  simulation  of  Experiments  1  and  2  are 
shown  in  Figure  10.  (Variations  were  made  to  the  O(e)  relation  directly;  the  a(5)  relation 
was  then  derived.)  The  strength  (peak  of  the  stress-strain  curve)  assigned  to  the  rod  was 
4.4  MPa,  30%  higher  than  the  measured  static  splitting  strength.  In  both  cases,  the 
strength  at  the  complete  fracture  location,  (0.46  cm  from  the  midpoint  in  Experiment  1  and 
0.1  cm  from  the  midpoint  in  Experiment  2)  was  degraded  to  3.5  MPa  (80%  of  the  strength 
in  the  rest  of  the  rod)  to  produce  complete  fracture  only  at  the  desired  location.  The  critical 
fracture  volume  per  unit  area  was  about  6  Jim.  This  is  slighdy  larger  than  one  might 
expect;  it  corresponds,  for  example,  to  the  volume  of  a  5-cm-diameter  crack  under  a  free- 
field  tension  of  about  3  MPa. 

The  strength  levels  were  chosen  to  obtain  a  good  prediction  of  the  peak  strains. 
(Higher  strength  gives  higher  peak  strains.)  The  indtial  softening  slope  of  the  stress-strain 
relation  was  chosen  to  match  the  rounded  peaks  in  the  strain  histories.  (A  steeper  slope 
gives  more  pointed  peaks,  a  less  steep  slope  gives  flatter  peaks.)  The  critical  fracture 
volume  per  unit  area  was  chosen  so  that  the  pulse  duration  of  the  computed  strains  would 
match  that  of  the  measured  strains,  and  to  prevent  multiple  fractures.  (A  smaller  critical 
value  results  in  a  shorter  duration  pulse  and  multiple  fractures.)  The  1  MPa  difference 
between  the  two  curves  out  to  the  ultimate  strain  was  needed  to  prevent  multiple  fractures. 
(A  smaller  difference  would  result  in  multiple  fractures.)  Although  the  curves  in  Figure  10 
were  chosen  strictly  to  produce  a  good  match  to  the  measured  strain  histories,  it  is  interest¬ 
ing  that  the  form  of  the  curves  is  similar  to  the  static  measurements  made  by  Wiliam,  et  al 
(1984)  on  a  lower  strength  concrete. 

Experiment  1 

The  comparisons  of  calculated  and  measured  strains  at  ±7.62  cm  in  Experiment  1  are 
shown  in  Figure  11.  (The  calculated  strains  for  Experiments  1  and  2  are  the  average  of  the 
four  cells  spanning  the  strain-gage  locations.)  As  the  comparisons  of  strain  histories  show, 
this  set  of  stress-strain  relations  produces  a  very  good  match  to  the  experimental  results. 
The  calculated  peak  strains  for  this  case  are  shown  in  Figure  12.  At  the  locations  of  the 
strain  gages,  the  calculated  inelastic  strain  (averaged  over  the  length  of  the  gage)  was  only 
about  10  microstrain.  Thus,  the  assumption  that  the  measured  peak  strains  were  elastic 
would  give  a  reasonable  approximation  to  the  strength.  However,  inelastic  strain  occurred 
nearly  everywhere  in  the  rod,  so  a  strictly  elastic  analysis  would  not  be  appropriate.  In 
particular,  damage  was  concentrated  at  two  locations  about  3  cm  from  the  fracture.  The 
magnitudes  of  these  peaks  imply  that  less  than  half  of  the  unrecovered  work  done  was 


Measured 

Calculated 


r 


Calculated 

Measured 


1  — i - i - *—  1  ■  i  —  J  —  <  .  ■■  i  »  - 

0.25  0.3  0.35  0.4  0.45  0.5 

TIME  (ms) 

pm  :  i 

(b)  X  =  7.66  cm 

JA-4451  -95B 

IW\ 

• 

Figure  1 1 .  Calculated  and  measured  strains  tor  Experiment  1 . 

1 

mm  •••-■■■! 

W 

B-18 

T 

ftji 

'I'ffllSS? 

SJk!K#C«2 

•SwSwBWW 

jgjggj 

wwvSSS? 

Fracture  Location 


JA-4451  -9oA 


Figure  1 2.  Peak  strains  calculated  tor  Experiment  1 . 
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energy  released  at  the  location  of  complete  fracture. 

Experiment  2 

To  simulate  Experiment  2,  the  same  stress-strain  relations  were  used,  but  the  weak¬ 
ened  cell  was  centered  at  the  midpoint  of  the  rod.  (In  the  experiment,  the  fracture  location 
was  only  0.10  cm  from  the  midpoint  of  the  rod.)  Two  simulations  were  performed,  one 
with  linear  unloading  to  the  origin,  and  one  with  elastic  unloading,  as  shown  in  Figure  10. 
The  results  are  compared  with  the  strains  measured  in  the  experiment  in  Figure  13.  The 
strain  measured  at  +7.62  cm  compares  very  well  with  the  strain  predicted  using  inelastic 
unloading.  On  the  other  hand,  the  strain  measured  at  -7.62  cm  has  a  residual  amplitude 
similar  to  the  prediction  with  elastic  unloading.  No  attempt  was  made  to  improve  the  com¬ 
parisons  by  changing  the  stress-strain  relations. 

Because  the  loading  and  fracture  location  in  this  experiment  were  essentially  sym¬ 
metric,  the  strains  at  ±7.62  cm  would  be  expected  to  be  the  same.  Because  the  measured 
strains  were  different  at  the  two  locations,  the  effects  of  inhomogeneities  not  included  in 
the  simulation  must  have  caused  this  asymmetric  response.  The  two  simulations  with 
different  unloading  paths  suggest  that  the  asymmetry  is  at  least  partly  associated  with  inho- 
mogeneiuies  in  unloading  characteristics. 


DISCUSSION 

The  first  point  of  discussion  is  the  fact  that  this  trial-and-error  exercise  certainly  does 
not  produce  a  unique  interpretation  of  these  experiments.  Within  the  assumptions  of  the 
model  however,  that  the  failure  is  associated  with  cracking  and  that  the  cracks  are 
separated  by  finite  dimension,  we  found  that  there  is  a  very  narrow  range  of  parameter 
choices  that  produce  agreement  with  the  experimental  observations.  That  is,  seemingly 
small  deviations  from  the  parameters  given  here  result  in  much  poorer  agreement  with  the 
experiments.  The  critical  test  of  the  parameter  set  lies  in  matching  the  strain  histories 
while  also  preventing  multiple  fractures. 

Secondly,  whereas  the  dynamic  experiments  produced  a  single  complete  fracture  with 
no  visible  secondary  damage,  the  calculations  predicted  some  inelastic  strain  to  occur  virtu¬ 
ally  throughout  the  specimen,  with  concentrations  of  inelastic  strain  within  a  few  centime¬ 
ters  of  the  locations  of  complete  fracture.  Simulations  in  which  this  distributed  damage  did 
not  occur  showed  poor  agreement  with  the  measured  strains.  This  result  suggests  that 
extensive  microcracking  in  the  concrete  may  have  occurred  in  the  experiments  in  the  region 
of  the  specimen  surrounding  the  complete  fracture  location.  In  fact,  a  posttest  microscopic 
inspection  made  by  Gran  and  Seaman  (1986)  of  the  interior  of  a  specimen  from  a  similar 
experiment  did  show  distributed  microcracking  that  is  qualitatively  consistent  with  the 
strain-softening  predictions  of  maximum  strain.  These  microcrack  observations  have  not 
yet  been  fully  quantified,  however,  so  a  direct  comparison  is  not  shown  here. 
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Finally,  we  noticed  that  the  calculated  stress  at  the  complete  fracture  location  (not 
shown)  drops  from  3.5  MPa  to  zero  in  about  10  (is.  This  observation  raises  the  question 
of  how  the  strain  history  at  the  measurement  locations  can  have  a  decay  time  of  about  50 
(is.  The  answer  is  that  the  inelastic  strain  between  the  fracture  location  and  the  strain  gage 
locations  drastically  alters  the  tensile  stress  wave  as  it  propagates.  In  fact,  for  matching  the 
strain  measurements,  the  description  of  the  material  behavior  in  this  region  has  a  far  greater 
effect  than  does  that  at  the  location  of  complete  fracture.  (A  calculation  in  which  the  cell 
at  the  fracture  location  is  given  the  same  strength  as  before  but  with  only  half  the  ultimate 
strain  produces  a  much  quicker  decay  of  stress  at  the  fracture  location  but  only  a  slightly 
quicker  decay  in  strain  at  the  strain  gage  locations.)  That  is,  the  calculations  are  not  as 
sensitive  to  the  details  of  the  softening  properties  at  the  location  of  complete  fracture, 
except  that  fracture  must  occur  there.  Thus,  the  properties  assumed  for  the  rest  of  the  rod, 
rather  than  those  assumed  for  the  fracture  location,  are  a  better  description  of  the  material 
in  general. 

SUMMARY  AND  CONCLUSIONS 

A  one-dimensional  strain-softening  model  was  used  in  wave-propagation  calculations 
to  interpret  results  of  dynamic  tension  experiments  on  concrete  rods.  The  model  is  based 
on  the  assumption  that  the  stress-strain  relation  is  not  a  property  of  a  material  point  (as  in 
continuum  theory),  but  is  an  average  property  of  a  finite  volume  of  material  containing  a 
developing  crack  or  failure  plane.  The  stress-strain  relation  thus  has  associated  with  it  a 
finite  dimension,  namely  an  effective  crack  separation  distance. 

Using  the  idealized  one-dimensional  strain-softening  stress-strain  relation,  two 
dynamic  unconfined  tension  experiments  were  simulated  several  times  with  trial  sets  of  the 
material  parameters  until  good  agreement  with  the  measured  strains  was  obtained  in  both 
cases.  The  two  experiments  were  simulated  with  the  same  set  of  stress-strain  relations, 
which  is  consistent  with  the  fact  that  the  experimental  results  were  very  similar  both  qualti- 
tatively  and  quantitatively.  In  addition  to  providing  an  estimate  of  the  dynamic  tensile  pro¬ 
perties  of  the  concrete,  these  calculations  suggest  that  tensile  damage  in  the  concrete  was 
distributed  over  several  centimeters.  Finally,  the  calculations  suggest  that  the  strain  history 
measured  a  few  centimeters  from  the  location  of  complete  fracture  is  primarily  a  function 
of  inelastic  wave  propagation  from  this  location  to  the  strain  gage  (through  a  region  of  dis¬ 
tributed  tensile  damage),  and  is  less  dependent  on  the  behavior  of  the  material  at  the  loca¬ 
tion  of  complete  fracture. 

The  stress-strain  relations  used  to  obtain  the  best  match  with  the  data  certainly  cannot 
be  considered  unique,  but  there  is  not  very  much  latitude  in  the  choice  of  the  material 
parameters  that  produce  a  good  match  with  the  measurements.  Based  on  these  analyses, 
the  unconfined  tensile  strength  of  the  concrete  at  a  strain  rate  of  10/s  to  20/s  is  about  4.4 
MPa,  nearly  30%  higher  than  the  static  splitting  tensile  strength  of  3.4  MPa.  The  stress 
versus  fracture  volume  relation  is  not  linear,  the  critical  fracture  volume  per  unit  area  is 


about  6  Jim,  and  the  effective  crack  spacing  is  about  0.635  cm.  The  critical  fracture 
volume  is  larger  and  the  effective  crack  spacing  smaller  than  would  be  derived  from  static 
experiments  and  analyses.  These  parameters  are  apparently  dependent  on  the  dynamics  of 
fracture. 
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APPENDIX  II.  -  NOTATION 

The  following  symbols  are  used  in  this  paper: 

E  *  elastic  modulus 

G  =  stress 

Gq  =  tensile  strength 

E  =  strain 

Eg  =  strain  at  peak  stress 

£ult  =  ultimate  strain 

lQ  =  initial  cell  dimension 

5  =  fracture  volume  per  unit  area 

8C  =  critical  fracture  volume  per  unit  area 

Gc  =  fracture  energy  per  unit  area 

G'  =  work  done  to  partially  fracture  a  cell 


APPENDIX  C 

A  SUBROUTINE  FOR  DETERMINING  STRESS  INTENSITY  VALUES, 
CRACK  OPENINGS,  AND  EFFECTIVE  MODULI  FOR 
MULTIPLE  CRACKS  IN  AN  ELASTIC  MATERIAL 

1.  INTRODUCTION 

Kachanov1  has  constructed  an  approximate  analysis  for  computing  stress  inten¬ 
sity  values  for  arrays  of  cracks  of  arbitrary  locations,  lengths,  and  orientations  in 
an  elastic  solid.  This  report  outlines  an  implementation  of  Kachanov’s  procedure 
into  a  computer  program  for  routine  evaluation  of  these  stress  fields. 

The  method  is  based  on  classical  solutions  for  the  stress  states  around  a  two- 
dimensional  flat  crack  in  an  infinite  elastic  body  under  external  tractions.  The 
procedure  provides  for  a  simultaneous  solution  for  the  stress  intensity  values  [Kr 
and  Ku)  at  each  end  of  each  crack  in  the  array.  With  these  K  values,  the  crack 
opening  shape  is  determined.  From  the  crack  shape,  the  average  crack  opening 
strain  in  the  body  can  be  found. 

This  Appendix  outlines  the  basic  analytical  procedure  developed  by  Kachanov 
and  then  describes  the  techniques  used  in  our  program.  Some  comparisons  with 
exact  solutions  are  given. 

2.  OUTLINE  OF  KACHANOV’S  METHOD 

Kachanov  considers  an  array  of  two-dimensio,  cracks  embedded  in  an  elastic 
material  under  a  remote  loading  (<7°°).  We  assume  that  the  cracks  are  of  arbitrary 
length,  location,  and  orientation.  The  length  of  the  rth  crack  is  given  by  Lr  and 
the  orientation  by  the  unit  normal  nr. 

To  quantify  these  cracks,  Kachanov2  forms  a  crack  density  tensor  A  with  the 
following  definition: 

C-] 


and  represents  the  normals  in  the  i,  j  coordinate  system  thus: 


nr  =  +  n^j 


nr  =  cos  6ri  +  sin  0rj  =  —  sin  <j>ri  +  cos  4>rj 


where  6r  is  the  angle  from  the  i  direction  to  the  normal  to  the  crack  plane  and 
4>r  is  the  angle  from  i  to  the  crack  plane.  Note  that  the  factor  nrnT  on  the  right 
in  Equation  (1)  is  a  dyad,  not  a  vector  or  dot  product.  Then  the  components  of 
the  crack  density  tensor  are 


A<i  =  £  iWM” 


From  its  definition,  we  can  see  that  A  is  a  symmetric  second  order  tensor  (as 
a  sum  of  such).  This  tensor  describes  the  lengths  and  orientations  of  the  cracks 
but  not  their  locations  in  the  material.  The  principal  axes  of  A  are  the  axes  of 
orthotropy.  In  more  detail,  A  is  given  by  the  following  equation: 


A  =  '52  L2r  [sin2  4>ri  i  +  cos2  <j>rj  j  -  sin  <f>r  cos  <j>T{i  j  +  j  i')] 


or  the  matrix  of  components 
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Er  Ll  sin2  <t>r 


-  Er  L2r  sin  <l>r  COS  <j>r 


Er  Sin  4>r  COS  4>r 
Er^COS2^ 


The  stress  acting  along  each  crack  is  a  superposition  of  the  external  stress  field 
on  the  stress  field  caused  by  the  opening  of  each  of  the  other  cracks.  To  proceed, 
Kachanov  considers  just  two  cracks  and  their  interactions.  He  evaluates  first  the 
interaction  stresses  cr"(£)  and  ajr(£):  the  stresses  at  position  £  along  the  sth  crack 
generated  by  uniform  normal  and  shear  tractions  of  unit  intensity  along  the  rth 
crack.  To  obtain  the  actual  stresses  along  the  sth  crack,  we  multiply  by  the  average 
of  the  normal  or  shear  stresses  on  the  Tth  crack:  <  Pr  >  and<  rT  >.  This  use  of 
the  average  stresses  on  the  cracks  is  the  essential  ingredient  of  Kachanov’s  method 
that  simplifies  the  solution  enough  to  make  it  practical. 


With  the  foregoing  definitions,  we  can  now  write  the  basic  equation  for  the 
normal  and  shear  stress  at  any  point  £  along  the  sth  crack  as  a  function  of  the 
external  normal  and  shear  loading  (Pf°  and  rf5)  and  the  interaction  stresses  from 


the  other  cracks. 


p.( 0  =  Pr  +  «.[£<*(«)  <  P,  >  +<C(f)  <  n  >]«. 


*■.(«)  =  rr  +  n.\E<{()  <Pr>  +°'M)  <  T,  >1  -  fj>,({)  -  P“ In,  (S) 


Here  P3°°  —  n,a°°  n,  and  rs°°  —  n,a°°  —  Pf°  n ,  are  the  normal  and  shear  tractions 
induced  along  the  sth  crack  length  by  the  remote  loading  in  the  absence  of  the 
cracks.  The  n ,  factors  are  normals  to  the  plane  of  the  sth  crack;  hence,  their 
presence  in  the  equations  provides  for  a  standard  angular  transformation  of  the 
applied  stresses  onto  the  sth  crack  plane.  The  quantities  cr"r  and  oT  are  the  stresses 
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generated  along  the  sth  crack  trace  by  uniform  tractions  (normal  and  shearing)  of 
unit  intensity  acting  on  the  Tth  crack.  These  stresses  are  listed  in  Section  5  below 
for  convenience.  The  <  Pr  >  and  <  fr  >  quantities  are  obtained  by  averaging  the 
normal  and  shear  stresses  acting  along  the  rtfl  crack. 

The  first  step  toward  solution  is  to  average  the  stresses  in  Equations  (7)  and 
(8)  along  the  crack  lengths. 


<  P .  >=  pr  +  E(A”"  <  Pr  >  +a™  <  TV  >) 
r^aj 


<  T,  >=  r4°°  +  5Z(Cr  <  Pr  >  +A"  <  Tr  >)  (10) 

»#r 

The  A  quantities  in  these  equations  combine  the  averaging  process  over  the 
length  of  the  sth  crack  and  the  tensor  transformation  for  orientation. 


Kr  =n,<  cr;r  >  ft. 


Kr  =n,<  a*  >  -A (12) 
where  k  refers  to  either  n  or  r. 

The  second  step  in  the  solution  is  to  compute  the  A  quantities  in  Equations 
(9)  and  (10).  The  simultaneous  equations  given  by  Equations  (9)  and  (10)  can  be 
written  in  matrix  form  as  follows: 


TP  =  F 


C-A 


where  T  is  a  matrix  based  on  the  A  values,  P  is  a  vector  containing  the  averaged 
quantities  <  P,  >  and  <  ra  >,  and  F  is  a  vector  of  the  external  tractions  Pa°°  and 

T?- 

Next,  the  matrix  T  is  inverted  to  solve  for  P.  With  the  average  values  of  the 
stresses  on  each  crack  known,  we  can  return  to  Equations  (7)  and  (8)  to  compute 
the  stresses  Pa(£)  and  ra(£)  along  the  crack  lengths. 

Stress  Intensity  Factors 

With  the  stresses  known  along  the  crack  lengths,  the  stress  intensities  are 
computed  from  the  following  formulas: 

1  rL  \L  ±  ^l1/2 

*'(±i)  =  ^r/Jr?f]  p[m  (14) 

and 

1  fL  [L±  £1X/2 

^  <15> 

where  L  is  half  the  crack  length. 

Crack  Opening  Displacement 

The  crack  opening  (normal  displacement)  is  computed  using  the  following 
shape  function  for  a  quadratically  distorted  ellipse,  obtained  by  multiplying  the 
expression  for  the  coordinate  of  an  ellipse  by  a  quadratic.  By  using  the  quadratic, 
Kachanov  is  allowing  the  crack  shape  to  distort  from  the  ellipse  that  would  be 
obtained  for  an  isolated  crack.  The  normal  crack  opening  (6n)  is 

K  =  1  +  “4  +  V1  -  e!/£J  <16) 

A  similar  formula  is  written  for  bT,  the  shearing  displacement.  Kachanov  gives 


expressions  for  S„,  an ,  and  /?„  in  terms  of  Kj  and  Kjj  at  both  ends  of  the  crack. 


Sn  =  2  <  P,  > 


KjW+Kjj-L) 

2  \/nL 


(17) 


„  ml)  -  MzB. 

"  4  <P,>  yfil  -  [K^L)  +  Kj{-L)\ 

Kj{L)  +  Kj(-L)  -2  <P,>  VjFZ 
P  4  <  P,  >  \^L  -  [Kj{L)  +  M-L)) 

For  the  shearing  displacement  bT  identical  expressions  are  used  for  ST,  aT,  and 
0T ,  except  that  Kj  is  replaced  by  Ku  and  <  Ps  >  is  replaced  by  <  t,  >.  In  all 
these  foregoing  equations  for  the  crack  opening,  <P,>  and  <t,>  are  the  values 
from  Equations  (9)  and  (10);  that  is,  the  average  normal  and  shearing  stresses 
on  the  crack  face.  Hence,  with  the  K  values  known ,  the  normal  and  shearing 
displacements  of  the  crack  surfaces  can  be  found. 

The  average  opening  displacement  of  the  cracks  is  given  by  integrating  the  bn 
value  over  the  area  of  the  crack  face: 


<&„>=  ^(l+ft,/4) 


(20) 


and  the  crack  volume  (Bn)  is 


Bn  =  ~^(1  +  A./4) 


(21) 


Similar  expressions  can  be  formed  for  the  average  and  total  shearing  distortions 
<  bT  >  and  Br. 


Effective  Elastic  Moduli 


The  derivation  of  the  effective  moduli  begins  with  the  determination  of  the 
appropriate  coordinate  axes  for  the  modulus  tensor.  Kachanov  notes  that  if  crack 
interactions  are  neglected,  the  effective  properties  are  always  orthotropic  (that 
is,  they  have  rectangular  symmetry)  for  any  crack  configurations.  Thus,  in  the 
approximation  of  noninteracting  cracks,  the  axes  of  orthotropy  provide  a  natural 
coordinate  system  for  effective  moduli  when  interactions  are  taken  into  account. 
Therefore,  Kachanov  outlines  the  following  steps  for  determining  the  effective 
moduli: 

Step  1:  Find  the  axes  of  orthotropy  (?i,  e2)  assuming  the  cracks  do  not  interact. 

Step  2:  Find  the  effective  properties  for  interacting  cracks  in  the  coordinate 
system  ?x,  e2. 

In  the  first  of  these  steps  the  axes  of  orthotropy  axe  determined  using  the  A  matrix 
given  in  Equation  (6).  The  angle  a  measured  counterclockwise  from  the  i  direction 
in  the  t,  j  external  coordinate  system  is  computed  as  follows  from  the  components 
of  the  A  matrix: 


tan  2a  = 


(22) 


We  begin  step  2  by  looking  for  matrix  of  effective  compliances: 


<  tij  >—  CijkiCTu  (23) 

where  <>  indicates  a  volume  average  and  the  stresses  are  the  remotely  applied 
ones.  Hereafter,  t,  j ,  k,  and  l  refer  to  the  orthotropic  axes  t\  and  e2. 

The  general  formula  for  the  derivation  of  the  elastic  coefficient  matrix  Ciju  is 


<  £  >=  C°  :  a°°  +  ~  52(nrbr  +  brnr)L, 


(24) 
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where  5  is  the  area  of  the  elementary  cell  containing  the  cracks,  C°  is  the 
elastic  compliance,  and  ct00  is  the  externally  applied  stress.  The  n  vectors  are 
normal  unit  vectors.  br  is  the  average  displacement  discontinuity  across  the  rth 
crack.  Referring  to  Equation  (20),  br  has  the  definition 


br  =<  bn  >  mn+  <  bT  >  mT 


and  mn  and  mr  are  unit  vectors  in  the  direction  normal  and  along  the  rth  crack. 


For  br  Kachanov  makes  the  following  approximation  to  Equation  (20): 


h  -  11T  T  7 
or  —  -—LrtT 


where  tr  is  the  average  traction  provided  by  Kachanov’s  method  and 


and  7  is  assumed  to  be  1.  Here  Kachanov  is  neglecting  the  fact  that  crack 
opening  depends  on  the  orientation  of  the  applied  stress.  Now  Equation  (24)  can 
be  written 


<  £  >=  C°  :  2“  +  —  £(»,<;  +  t,n ,)L\ 


To  construct  the  fit  quantities  in  the  preceding  equation,  we  let 


tr  —  <rlel  +  tr2e2 


(29) 


Now  we  rewrite  Equation  (2)  as 

nT  =  ?!  +  n^e 2  (30) 

Then  we  combine  this  definition  of  nr  with  that  for  tT  and  form  the  products 
in  Equation  (28). 

ntr  =  tri(n{rl)eiei  +  n^eye.)  +  ^(ri^eyey  +  n^eyey)  (31) 

With  reversed  symbols,  the  equation  is 

trn  =  ?ri(^1)e,e,-  +  n(2)e,ey)  +  ?r2(n(1)eJe,-  +  n^2)e,e;)  (32) 

The  only  difference  between  Equations  (31)  and  (32)  is  the  reversal  of  the  vector 
pairs  e,- e y  in  the  central  two  terms.  With  Equations  (31)  and  (32)  in  Equation  (28), 
we  can  see  the  form  the  compliance  equation  must  take.  The  stress  quantities  crkl 
in  Equation  (23)  are  applied  one  at  a  time,  so  the  tr  quantities  always  represent 
just  one  applied  stress  component.  The  vectors  e,  and  ey  then  correspond  to  the 
directions  of  the  strain  e,y . 

Let  us  now  examine  the  nt  expressions  in  preparation  for  programming  them 
into  the  code.  First  examine  the  nr  terms  in  Equation  (30). 


nr  —  cos(0r  —  a)e i  4-  sin(0r  —  a)e2 


(33) 


The  £rl  and  tk2  terms  are  derived  in  a  three-step  procedure: 

1.  Apply  a  unit  external  stress  in  the  ei-e2  orthotropic  plane,  and  transform 
this  stress  to  the  fixed  i-j  plane  to  obtain  the  components  crZI,  ayv,  and  axy. 

2.  Determine  the  normal  and  shear  stresses  on  each  plane  from  this  external 
loading. 

3.  From  the  Pj50  and  t”  arrays  and  the  T  matrix,  derive  the  modified  normal 
and  shear  stresses  (<  P,  >  and  <  t,  >)  on  each  crack,  accounting  for  the 
presence  of  the  other  cracks. 

Now  the  argument  in  the  sum  in  Equation  (28)  can  be  written  out.  The  expression 
is 

nrtr  +  trnr  =  2tri  cos(0r  —  a)eifii  +  (£ri  sin(0r  —  a)  +  £r2  cos(0r  —  a))e2<?i 

+  (frisin(0r  —  a)  +  fr2cos(0r  —  a))eie2  +  2£r2sin(0r  -  a)e2e2(34) 

Hence,  with  each  loading  o-jtf,  we  can  derive  4  e  terms  (en,  e22,  and  ei2  = 
e2i).  The  applied  external  stresses  are  an,  a22,  and  a12  (=  cr2i).  Note  that  in  our 
transformations  we  recognize  that  <7i2  and  cr2i  are  always  applied  together. 

Kachanov  begins  the  computation  of  the  compliances  with  the  application 
of  a  ‘trial’  stress  He  states  that  the  <  en  >  response  to  this  stress  provides 

the  effective  Young’s  modulus  (Pi)  in  the  directions,  or  Cun. 

a°°  =  £7^?!? i  (35) 

and  we  can  take  ct”  35  unity.  Then 


Here  tr\  is  the  component  of  tr  along  ef  calculated  under  conditions  of  our 
loading  off.  Thus 


E°  l+££rL?<trl> 


Similar  computations  will  yield  62222- 


The  computation  for  C2211  proceeds  as  follows.  Let  the  applied  stress  be  the 
same  as  in  Equation  (35)  and  take  <7“  as  unity.  Then 


<  *22  >~  -TT  +  cH2'C'^r2«r2)  =  C2211  =  ~  JT 

£j0  L£j0O  r  Ej\ 


where  tr 2  is  the  e2  component  of  tr  calculated  under  conditions  of  our  loading 


Next  we  compute  Cm 2,  the  <  612  >  response  to  off. 


The  trl  and  tr2  values  are  components  of  tr  as  calculated  under  loading  off  = 


The  constants  written  above  are  the  only  non-zero  ones  in  orthotropic  material. 
The  rest  of  the  constants  developed  below  characterize  deviation  from  orthotropy. 
For  these  constants,  let  us  start  with  Cmi,  the  <  ei2  >  response  to  off  =  1. 


<£,j>=2£??  «  =  C, 


The  tr  components  are  calculated  under  the  loading  of  crjf 


1.  We  can 


compute  C1222  similarly. 

Moduli  for  Noninteracting  Cracks. 

These  moduli  are  for  comparison,  to  estimate  the  effect  of  crack  interactions 
on  the  overall  moduli.  They  are  given  by  the  same  formulas  as  the  corresponding 
moduli  for  the  interacting  cracks,  the  only  difference  being  that  the  tractions  on 
cracks  ( tr )  are  taken  as  directly  induced  by  the  remote  loading  a°°;  that  is, 


tr  =  nr(l° 


(41) 


The  factor  on  the  right  is  a  dot  product  between  the  vector  nr  and  the  stress 
tensor.  Thus,  in  the  formula  for  Cun,  for  example, 


t,  =  n,  (a^ejfi)  =  nl'1?. 


(42) 


and  therefore, 


tri  =  and  tr 2  =  0 
and  hence  the  normalized  modulus  Ei/ E0  is 


Ei  =  1 

E0  1  +  %ZrL*(nl1))* 


(43) 


(44) 


and  similarly  for  C22 22- 

Other  compliances,  such  as  C2211,  are  unaffected  because  tri  is  zero  under 
loading.  That  is, 


For  C 1212  the  expression  is 


tr  -  nTof2{e ie2  +  c2ei) 
=  n^e2  +  n(2)ei 
=  tr2?2  +  trl  Cl 


Then  C1212  is  computed  from 


~  2ST  +  ^  [(n^)2  +  (n-2)^]  Lr 


—  4  I  7T  r  2  v  7 

The  foregoing  procedure  is  implemented  in  the  computer  program  to  be  de¬ 
scribed  next. 


3.  DESCRIPTION  OF  THE  KCRACK  PROGRAM 

The  calculations  in  the  program  begin  with  the  insertion  of  input  quantities 
to  describe  the  problem.  The  material  constants  E  (Young’s  modulus)  and  u 
(Poisson’s  ratio)  are  inserted  first.  Then  for  each  crack  the  crack  length  ( L ),  crack 
center  (X0,  Y0),  and  angle  (5)  of  the  normal  to  the  crack  are  inserted.  The  external 
field  is  described  by  Sxx,  Svv,  and  Sxv. 

The  normal  and  shear  stresses  on  the  cracks  caused  by  the  external  loadings 
are  determined  by  a  standard  tensor  rotation  in  the  subroutine  CRFORC: 


P,00  =  Sxx  cos2  9,  +  Syv  sin2  9,  +  2 Szy  sin  9 ,  cos  6, 


C-13 


r4°°  =  ~{SXX  -  Svv)  sin  Oa  cos  0a  +  Siy[cos2  9,  -  sin2  9, 


(49) 


These  quantities  —  P4°°  and  r4°°  —  are  contained  in  the  FORCE  array.  A 
similar  tensor  transformation  provides  for  evaluating  the  A  quantities: 

A.?  =  [««*(#,  -  *,)  /t_  (50) 

+2sin(04  -  9r)cos(93  -  6r)  f  sr{£)dt  +  sin2(04  -  9r)  [  <7yy4r(£)d£ 

The  oxxar  quantity  is  the  stress  in  the  X  direction  (of  the  rlh  crack)  from  the 
standard  stress  fields  for  normal  stress  on  the  rth  crack  acting  at  the  location  £ 
on  the  s‘h  crack.  In  constructing  this  stress,  we  use  the  standard  fields  given  by 
Kachanov  and  listed  in  Section  5,  the  distance  between  the  center  of  the  rth  crack 
and  the  point  £  along  the  sth  crack,  and  the  relative  orientations  ( 9r  and  9„)  of 
tne  two  cracks.  The  same  expression  is  used  for  A4IT  by  replacing  the  on  quantities 
with  oT . 

The  corresponding  expression  for  A™  is 

Kr  =  JT  Sin ^  COS(d>  ~  (JL  GXZ,rU)dC  +  fL  (0^)  (51) 

^cos2 (04  -  er)  -sin2(04  -  9r)^j  fL  cr"„4r(£)d£j 

A  similar  expression  is  used  for  A4^,  with  on  replaced  by  oT . 

The  integrals  in  the  foregoing  expressions  are  evaluated  by  Simpson’s  rule  in 
the  program  in  the  CRSTRS  subroutine.  The  integration  was  tested  using  10,  20, 


C-14 


40,  and  100  intervals  for  the  case  of  two  parallel  end-to-end  cracks  with  a  spacing 
of  0.01  L,  where  L  is  half  the  crack  length.  For  40  intervals  or  more,  the  results 
were  accurate  to  four  significant  figures. 

With  the  A  values  known,  we  construct  the  T  matrix  with  the  following  defi¬ 
nitions  for  the  components: 

Tkm  =  1  for  k  —  m  (52) 

=  0  for  fc  =  m  +  1  or  m  —  1 

=— A”"  for  r  —  kj 2  —  1/2  for  k  odd,  and  s  =  m/2  —  1/2  for  m  odd 
=  — A™  for  r  =  kj 2  —  1/2  for  k  odd,  and  s  =  m/2  for  m  even 
=  —  A"/  forr  =  kj 2  for  k  even,  and  s  =  m/2  —  1/2  for  m  odd 
=  —A”  for  r  =  k/2  for  k  even,  and  s  =  m/2  for  m  even 

This  T  matrix  (called  AMAT  in  the  program)  is  inverted  in  CRSTRS  to  obtain 
the  P  force  vector.  The  P  vector  contains  the  average  stress  (normal  and  shear) 
on  each  crack  (Equations  (9)  and  (10)). 

With  the  P  force  vector  known,  we  can  return  to  Equations  (7)  and  (8)  to 
compute  the  local  stresses  and  then  evaluate  the  stress  intensities  from  Equations 

(14)  and  (15).  These  local  stresses,  Pr{Z)  and  rr(£),  are  computed  in  the  subroutine 
CRKl  with  the  help  of  TRACTN.  Because  the  integrals  in  Equations  (14)  and 

(15)  are  singular  at  both  ends  of  their  range  of  integration,  we  used  Chebyshev- 
Gauss  quadrature  to  evaluate  them.  Equation  (14)  was  rewritten  as  follows  for 
the  numerical  integration: 


K,(L) 


1 

y/nL 

L 

y/nL 


a 

L 


L  +  Z 


1/2 


-lVL-Z  j 
(1  +  Z/L)P(t) 


yji  -  e/v 


P(Z)dZ 
d(Z/L ) 


(53) 


Then  the  numerical  expression  is 


K^L)Slt^~fl(1  +  Xi)P(Xi)  (54) 

m  I 

where  x,  =  cos[(2i  —  l)7r/(2m)],  and  m  is  the  number  of  intervals  used  in  the 
integration.  The  integration  for  K  was  tried  with  various  values  of  m.  Results 
that  were  accurate  to  four  figures  were  obtained  with  m  =  40. 

4.  EXACT  SOLUTIONS  FOR  STRESS  INTENSITY  FACTORS 

To  evaluate  the  accuracy  of  the  preceding  method,  it  is  useful  to  have  exact 
solutions  for  some  special  cases.  Such  solutions  are  available  for  the  case  of  two 
colinear  cracks.  In  addition,  Kachanov  was  able  to  obtain  an  analytical  solution 
based  on  his  approximate  method.  These  two  types  of  analytical  solutions  are 
presented,  and  some  representative  values  of  Ki  are  given. 

In  all  cases,  the  analysis  is  for  a  pair  of  cracks  extending  along  the  x  axis  from 
—  1  to  —  k  and  from  k  to  I.  The  elastic  body  is  loaded  by  a  uniform  tension  in  the 
Y  direction.  The  resulting  FT/  and  Ku  values  are  given  in  Table  C-I.  The  K  values 
listed  in  the  columns  labeled  “Procedure”  were  obtained  by  the  numerical  method 
discussed  in  this  report;  those  in  the  column  labeled  “Approx.”  by  an  analytical 
solution  to  the  method  reported  here,  and  given  by  Kachanov  in  his  paper;  and 
those  in  the  column  labeled  “Exact"  by  the  exact  analytical  solution.  The  results 
in  the  table  show  that  the  numerical  procedure  gives  an  acceptable  approximation 
to  the  exact  results  throughout  the  range  of  interest.  Also,  the  numerical  method 
is  able  to  reproduce  the  values  of  the  analytical  solution  to  this  procedure  within 
0.1%  for  crack  separations  down  to  1%  of  the  crack  length.  The  surprising  result 
of  Kachanov’s  procedure  is  that  it  is  so  accurate  for  very  close  spacings  of  the 
cracks. 

Exact  Solution  of  the  Approximation 

For  two  cracks  along  a  line,  an  analytical  solution  was  obtained  by  Kachanov 
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TABLE  C-I 

COMPARISON  OF  STRESS  INTENSITY  VALUES 


Kj(k) 

Errors,  % 

Procedure 

Approx. 

Exact 

Proc. 

Approx. 

0.200 

1.112018 

1.112018 

1.112470 

0.040 

0.040 

0.100 

1.250944 

1.250944 

1.255122 

0.333 

0.333 

0.070 

1.346863 

1.346862 

1.356894 

0.739 

0.739 

0.050 

1.452423 

1.452421 

1.472882 

1.389 

1.389 

0.020 

1.808529 

1.808426 

1.904569 

5.04 

5.05 

0.010 

2.134548 

2.133550 

2.371571 

9.99 

10.04 

0.007 

2  317831 

2.315134 

2.671641 

13.24 

13.34 

0.005 

2.499528 

2.493311 

2.999207 

16.67 

16.87 

0.002 

3.041578 

3.002006 

4.164502 

26.96 

27.91 

0.001 

3.517141 

3.399714 

5.394657 

34.80 

36.98 

■i 

*>(i) 

Errors,  % 

Procedure 

Approx. 

Exact 

Proc. 

Approx. 

0.200 

1.051580 

1.051580 

1.051682 

0.00965 

0.00963 

0.100 

1.085775 

1.085775 

1.086335 

0.0516 

0.0516 

0.070 

1.102807 

1.102807 

1.103874 

0.0963 

0.0967 

0.050 

1.118012 

1.118011 

1.119791 

0.1588 

0.1590 

0.020 

1.153792 

1.153772 

1.158939 

0.444 

0.446 

0.010 

1.174985 

1.174831 

1.184110 

0.778 

0.784 

0.007 

1.184115 

1.183738 

1.195629 

0.963 

0.995 

0.005 

1.191838 

1.191042 

1.205669 

1.147 

1.213 

0.002 

1.210363 

1.206282 

1.229385 

1.547 

1.879 

0.001 

1.224751 

1.214200 

1.244326 

1.573 

2.421 

Note:  The  “Procedure”  values  were  generated  using  40  intervals. 


for  the  Kj  values.  The  cracks  extend  from  —1  to  —  k  and  from  A:  to  1.  So  the  K 
values  have  the  following  expressions: 


r 


a 


IH 


=  1  +  (1;  _  A)t(1  -  fc)  [2B(m)  “  *1*  +  W™)  -  f (1  -  *)]  (55) 


K,(k)/K,.  =  1  +  [— 2£(”>)  +  (*  +  1  )*(»>)  -  1(1  -  *)]  (56) 


where  Ki0  =  o\J 7r(l  —  k)/ 2  is  the  stress  intensity  factor  for  an  isolated  crack, 

A  =  ^/2(l  +  fc)/(l  +  y/k)  —  1  is  the  transmission  factor, 
m  =  \/l  ~  k2  is  the  argument  of  the  elliptic  integrals, 
and  K  and  E  are  complete  elliptic  integrals  of  the  first  and  second  kind. 


These  expressions  for  Kj  were  evaluated  for  a  range  of  k  values  and  the  results 
are  listed  in  Table  C-I. 


Analytical  Solution 


Page  46  of  Sneddon  and  Lowengrub4  gives  formulas  for  Kj  and  Ku  values  for 
a  pair  of  adjacent  horizontal  cracks  under  either  normal  ( a )  or  shear  ( r )  stress 
applied  remotely.  The  cracks  extend  from  x  =  —  b  to  —k  and  from  k  to  b.  The 
expressions  for  the  K  values  at  the  near  ends  are 


K,w  -  0(1, Ik)  jpZTpyfi 


K„(k)  -  r(*/fc)  (JJ 


and  the  K  values  at  the  remote  ends  are 


TOE 


NVi 


*,(*)-,(»*)*> zirnm 

k 


Kn(i)  =  r^ynlzmLm 

K 


In  these  expressions,  K(m)  and  E(m)  are  complete  elliptic  integrals  of  the  first 
and  second  kind,  and  m  is  the  modulus: 


m  —  1  —  k2/b2  (61) 

For  comparisons  with  the  other  computations,  it  is  useful  to  display  these  stress 
intensity  factors  as  ratios  of  the  present  values  divided  by  the  K  values  for  isolated 
cracks.  These  isolated  K  values  are 

KIo  =  oJn{b-k)l  2  (62) 


and  a  similar  one  for  Kji0. 

These  expressions  for  Kj  were  evaluated  at  several  k  values  and  the  results  are 
listed  in  Table  C-I. 

5.  STANDARD  STRESS  FIELDS 

Kachanov3  has  provided  the  standard  stress  fields  (in  an  elastic  body)  resulting 
from  uniform  tractions  applied  along  the  face  of  a  crack.  For  these  equations,  the 
crack  lies  along  the  X-axis,  with  its  center  at  the  origin  of  coordinates.  The  uniform 
tractions  axe  p  and  r.  The  stresses  at  any  point  in  the  X-Y  plane  are: 

=  p{h  -  8Y2/4  +  8Y*I6) 

<7y  y  —  p(1 2  +  AY2 1\  —  8Y4Iq)  (63) 


WM®, 


•»  v  ' 
■>>>>> 


axv  =  2p{-Yh  +  XYh  +  AYzh-AXYih) 

Similarly,  for  the  shearing  stress: 

axx  =  2T{3Yh-3XYh-4Yzh  +  4XY3h) 

avv  =  2T{-YI3  +  XYh  +  4Y3h-4XYzh)  (64) 

oxv  =  t{I2-8Y2U+8Y*I6) 

The  I  factors  in  these  equations  are  defined  in  terms  of  the  locations  of  the 
point  of  interest  and  the  crack  length  as  follows: 


y/6(y/a  +  yFf  +  y/6) 


VS(V^  +  y/l  +  y/S) 

y/arjS3/2 

2i 

y/cFjS3/2 

L3  3y/arj{y/a  +  y/l)2{y/l  -  y/a)  +  <5(73/2  -  <*3  2) 


L 2  (a3'2  +  ~i2,2)6  +  3y/crj(y/a  + 


The  Greek  letters  are  defined  for  these  equations  as  follows: 

a  =  ( X-L)2  +  Y 2 
(3  =  2  {X2  +  Y2-L2) 

7  =  ( X  +  LY  +  Y 2 

6  =  p  +  2y/cFj 

These  are  the  standard  stress  fields  used  earlier  in  this  appendix. 


w*;< 


v 


v;v; 


6.  LISTING  OF  THE  KCRACK  COMPUTER  PROGRAM 


The  KCRACK  program  consists  of  a  main  program  plus  several  subroutines. 
The  following  listing  starts  with  KCRACK.  Next  is  the  COMMON,  a  listing  of  the 
shared  variables  for  use  in  the  main  program  and  all  the  subroutines  except  MINV, 
SIMQ,  and  TRACTN.  Following  the  COMMON  are  the  subroutines:  CRALF, 
CRCOMP,  CRFORC,  CRKl,  CRK2,  CRNINT,  CRSHAP,  CRSTIF,  CRSTRS, 
MINV,  SIMQ,  and  TRACTN.  The  purpose  for  each  subroutine  is  described  briefly 
in  the  subroutine  listing. 
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PROGRAM  KCRACK 

VERSION  OF  JANUARY  1988 
CHANGED  TO  HAVE  MANY  SUBROUTINES. 

PROGRAM  TO  COMPUTE  ELASTIC  CRACK  INTERACTION,  WRITTEN  BY 
LYNN  SEAMAN,  SRI,  JULY  1987  BASED  ON  DERIVATIONS 
FROM  MARK  KACHANOV,  TUFTS. 

X  AND  Y  ARE  COORDINATES  OF  -CENTROID  OF  CRACK, 

THCR  IS  THE  CCW  ANGLE  TO  THE  CRACK  NORMAL  IN  DEGREES. 
STRESSES  ARE  POSITIVE  IN  TENSION. 


INCLUDE  ' $DISK3 : [SEAMAN . CRACK] KCKCOM . FOR' 


CHARACTER* 10  Al , A2 , A3 , A4 , TITLE ( 8 ) 
DIMENSION  JBUG(IO) 

1-  KCRACK  3-  CRCOMP  5-  CRK1 

2-  CRALF  4-  CRFORC  6-  CRK2 


7-  CRNINT  9-  CRSTIF 

8-  CRSHAP  10-  CRSTRS 


VALUES  OF  JBUG ( 5 )  AND  JBUG(10)  OVER  2  GIVE  PRINTS  FROM  -TRACTN 


INPUT 


;  INPUT  ELASTIC  PROPERTIES 

;  EMOD  IS  YOUNG'S  MODULUS,  POISSON  IS  POISSON'S  RATIO 
:  ANINT  IS  THE  NUMBER  OF  INTEGRATION  INTERVALS  USED. 

PI  =  3.14159265358979323846 
WRITE  (6,1000) 

1000  FORMAT  (30X, ' S  RI  KCRACK'/'  CRACK  STRESS  INTENSITY 

1  'FACTORS,  CRACK  SHAPES,  AND  EFFECTIVE  STIFFNESSES'/) 

READ  (5,1006)  TITLE 
WRITE  (6,1006)  TITLE 

1006  FORMAT  (8A10) 

READ  (5,1001)  Al, EMOD, A2, POISSON, A3, ANINT, A4, JBUG 
WRITE  (6,1002)  Al, EMOD, A2, POISSON, A3, ANINT, A4, JBUG 
I BUG  *  JBUG ( 1 ) 

NINTRV  -  2 * INT (0.5 *ANINT  +  0.5) 

IF  (NINTRV  .GT.  100)  NINTRV  =  100 
IF  (NINTRV  .LT.  2)  NINTRV  =  2 
WRITE  (6,1010)  NINTRV 

1010  FORMAT  ('  ***  KCRACK  10,  POSSIBLE  RESETTING:  NINTRV  =  ',15) 

1001  FORMAT  (3 (A10,E10 .3)  ,  A10, 1011) 

1002  FORMAT  (3 (A10, 1PE10 .3) , A10, 1011) 

:  INPUT  EXTERNAL  STRESS  FIELD 

READ  (5,1003)  Al, NTYPE, A2 , SXX, A3, SYY, A4 , SXY 
WRITE  (6,1004)  Al, NTYPE, A2 , SXX, A3 , SYY, A4 , SXY 

1003  FORMAT  (A10 , 110 , 3 (A10 , E10 . 3) ) 

1004  FORMAT  (A10, 110, 3 (A10, 1PE10 . 3) ) 

:  INPUT  MICROCRACK  LENGTHS,  POSITIONS,  AND  ORIENTATIONS 
]  NCRACK  IS  THE  NUMBER  OF  CRACKS 

READ  (5,1003)  Al, NCRACK, A2 , SAREA 
WRITE  (6,1003)  Al, NCRACK, A2, SAREA 
DO  100  NC=1, NCRACK 

:  XCR  AND  YCR  ARE  CENTROID  OR  CRACK,  THCR  IS  ANGLE  TO  NORMAL 

READ  (5,1005)  Al, ELCRAK (NC) , A2 , XCR (NC) , A3 , YCR (NC) , A4 , THCR (NC) 
WRITE  (6,1007)  Al, ELCRAK (NC) ,A2,XCR(NC) , A3, YCR(NC) ,A4,THCR(NC) 

1005  FORMAT  (4 (A10, E10 . 3) ) 

1007  FORMAT  ( 4 ( A10 , 1PE10 . 3) ) 


INITIALIZATION  OF  XCRACK  AND  THETA 


THETA (NC)  =  PI/180 . *THCR(NC) 


Si 


S 


■a 


a 


9 


£ 


vl 


a 


c 


SINTH(NC)  =  SIN (THETA (NC) ) 

COSTH(NC)  =  COS (THETA (NC) ) 

ELSIN  =  ELCRAK(NC) *SINTH(NC) 

ELCOS  =  ELCRAK (NC) *COSTH (NC) 

XCRACK (NC,  1 )  =  XCR(NC)  -  ELSIN 
XCRACK (NC, 2 )  =  XCR(NC) 

XCRACK (NC, 3)  =  XCR(NC)  +  ELSIN 
YCRACK  (NC,  1 )  =  YCR(NC)  +  ELCOS 
YCRACK (NC, 2 )  =  YCR(NC) 

YCRACK (NC, 3)  =  YCR(NC)  -  ELCOS 

IF  (IBUG  .GT.  0)  WRITE  (6,  1098)  (XCRACK  (NC,  I)  ,  1  =  1,  3)  ,  XCR  (NC)  , 

1  (YCRACK (NC, I),I=1,3),YCR (NC) 

1098  FORMAT  ('  KCRACK  98  X=' , 1P3E15 . 8 , '  XCR=' , E15 . 8 /10X, ' Y=' , 3E15 . 8 , 
1  '  YCR=' , E15 . 8 ) 

100  CONTINUE 

C  ***  COMPUTE  THE  STRESSES  DUE  TO  THE  EXTERNAL  LOADING 
CALL  CRFORC ( JBUG ( 4 ) ) 


PART  1:  COMPUTATION  OF  THE  AVERAGE  STRESSES  ON  EACH  MICROCRACK 


CALL  CRSTRS ( JBUG (10)) 


PART  2:  DETERMINATION  OF  THE  -K-  STRESS  INTENSITY  FACTORS 


C  ***  DETERMINE  THE  TRACTIONS  ON  THE  CRACK  FACES 
CALL  CRK1  (JBUG (5)) 

C  ***  COMPUTE  THE  K  VALUES 
DO  700  NC  =  1 , NCRACK 
CALL  CRK2 (NC, JBUG (6) ) 


PART  3:  COMPUTATION  OF  THE  CRACK  SHAPES 


CALL  CRSHAP (NC, JBUG ( 8 ) ) 
700  CONTINUE 


PART  4:  COMPUTATION  OF  THE  COMPLIANCE  AND  STIFFNESS  MATRICES 


***  ALPHA  MATRIX  FOR  THE  ORTHOTROPIC  DIRECTIONS  FOR  NON-INTERACTING 
CRACKS . 


CALL  CRALF  ( JBUG (2) , SINALF, COSALF) 

C  ***  CONSTRUCT  SUM  (NT  +  TN)  LA2  FOR  Sll,  S22,  OR  S12 
CALL  CRSTIF  (JBUG ( 9 ), SINALF, COSALF) 

C  ***  CONSTRUCT  THE  COMPLIANCE  MATRIX 
CALL  CRCOMP ( JBUG ( 3 ) ) 

C  ***  CONSTRUCT  THE  STIFFNESS  MATRIX  FOR  INTERACTING  CRACKS 
CALL  MINV(AMAT, 6, DETERM, LLL,MMM) 

WRITE  (6,1908)  (AMAT(I) , 1=1, 36) 

1908  FORMAT  ('  STIFFNESS  MATRIX  =  ',1P6E12.5,'  FROM  KCRACK  908' 
1  /20X, 6E12 . 5/20X, 6E12 .5/20X, 6E12 .5/20X, 6E12 .5/20X,  6E12 .5) 


PART  5:  COMPUTE  COMPLIANCE  MATRIX  FOR  NON-INTERACTING  CRACKS 


CALL  CRNINT  (JBUG (7 ), SINALF , COSALF ) 

C  ***  CONSTRUCT  THE  STIFFNESS  MATRIX  FOR  NON-INTERACTING  CRACKS 
CALL  MINV(AMAT, 6, DETERM, LLL,MMM) 

WRITE  (6,1928)  (AMAT (I) , 1=1, 36) 

1928  FORMAT  ('  STIFFNESS  MATRIX  =  ',1P6E12.5,'  FROM  KCRACK  928' 

1  /20X, 6E12 . 5/20X, 6E12 .5/20X, 6E12 .5/20X, 6E12 . 5/2  OX, 6E12 .5) 

STOP  '  NORMAL  END  OF  KCRACK' 


IMPLICIT  REAL *8 (A-H, 0-Z) 

COMMONS  FOR  THE  -KCRACK-  CODE 
COMMON  KROW, NCRACK, NINTRV, SXX, SYY, SXY, EMOD, POISSON, SAREA 
COMMON  AK1 (10, 2) , AK2  (10, 2) 

COMMON  ELCRAK ( 10 )  ,XCR(10) ,YCR(10) ,XCRACK(10, 3) ,YCRACK(10, 3) 
COMMON  COSTH(IO) , SINTH(IO) ,THCR(10) , THETA (10) 

COMMON  SP (10, 101) , ST (10, 101) , SUMN(lOl) , SUMT (101) , XID (101) 
COMMON  SNN(10, 10, 101) , STN (10, 10, 101) , SNT (10, 10,101), 

STT (10, 10, 101) 

COMMON  ALPHA (2, 2) 

COMMON  AMAT (400) ,BMAT(400) ,FORCE(20)  ,PL(20)  ,TL(20,20> 

COMMON  SUM (3,3), COMPL (6,6), COMPLO (6,6) 

COMMON  SHAPES  (10, 2)  ,  SHAPEAdO,  2)  ,  SHAPEB  (10,  2) 

COMMON  LLL ( 6) , MMM (6) , PI 


c 

c 

c 


SUBROUTINE  CRALF  ( IBUG, SINALF, COSALF) 

COMPUTATION  OF  THE  ALPHA  MATRIX  TO  DETERMINE  THE  ORTHOTROPIC 
DIRECTIONS  FOR  NON-INTERACTING  CRACKS. 

INCLUDE  ' $DISK3 : [ SEAMAN . CRACK] KCKCOM . FOR' 

CRSUM  =  0 . 

DO  750  NC=1 , NCRACK 
EL2  =  ELCRAK(NC) **2 
CRSUM  =  CRSUM+EL2 
ALPHA (1,1)  =  ALPHA (1,1)  + 

ALPHA (1,2)  =  ALPHA (1,2)  + 

ALPHA (2,2)  =  ALPHA (2,2)  + 

750  CONTINUE 

ALPHA (2,1)  =  ALPHA(1,2) 

IF  (IBUG  .GT.  0)  WRITE 
1  ALPHA (2, 2) 

1754  FORMAT  (/'  ALF-754  ALPHA  =' , 1P2E12 . 5/ 12X, 2E12 . 5 ) 

DALPHA  =  ALPHA (1, 1) -ALPHA(2, 2) 

IF  (ABS'DALPHA)  .LT.  CRSUM* 1 . E-30 )  DALPHA  =  CRSUM*l.E-30 
TAN2ALF  =  -2 . *ALPHA ( 1 , 2 ) /DALPHA 
ALF  =  0 . 5*ATAN (TAN2ALF) 

ALFDEG  =  0 . 5*ATAN (TAN2ALF) *180/PI 
:  COS2ALF  =  0 . 5* (ALPHA(1, 1) -ALPHA(2, 2) ) /SQRT (ALPHAfl, 2) **2+ 

:  1  0 .25* (ALPHA(2, 2) -ALPHA (1, 1) ) **2) 

SINALF  =  SIN (ALF) 

COSALF  =  COS (ALF) 

]  -ALF-  IS  THE  ANGLE  OF  THE  ORTHOTROPIC  AXES  FROM  X,  Y,  (+  CCW) 

IF  (IBUG  .GT.  0)  WRITE  (6,1759)  ALFDEG, TAN2ALF, SINALF, COSALF 
1759  FORMAT  (/'  ALF-759  ALF (DEG) =' , 1PE12 . 5 , '  TAN2ALF=' , E12 . 5, 


EL2*SINTH (NC) **2 
EL2*SINTH (NC) *COSTH (NC) 
EL2  *COSTH (NC) **2 


(6,  1754)  ALPHA (1, 1) , ALPHA (1,2)  ,  ALPHA (2,  1) 


SINALF, COSALF='  ,2E12. 5) 


END 


SUBROUTINE  CRAMAT (AMAT, COMPL, POISSON) 
CONSTRUCT  THE  6X6  COMPLIANCE  MATRIX 


DIMENSION  AMAT (400) , COMP L ( 6 ,  6 ) 


DO  908  1=1,36 


908  AMAT ( I )  =  0. 

AMAT (1)  =  COMPL (1, 1) 

AMAT (2)  =  COMPL (1,2) 

AMAT (3)  =  -POISSON 
AMAT ( 4 )  =  2 . *COMPL (1,3) 

AMAT (7)  =  COMPL (2, 1) 

AMAT (8)  =  COMPL (2, 2) 

AMAT (9)  =  -POISSON 
AMAT (10)  =  2 . *COMPL (2,3) 

AMAT (13)  =  -POISSON 
AMAT (14)  =  -POISSON 
AMAT (15)  =  1. 

AMAT (19)  =  COMPL (3,1) 

AMAT (20)  =  COMPL (3,2) 

AMAT (22)  =  2 . *COMPL (3,3) 

AMAT (29)  =  2 . * (1 . +POISSON) 

AMAT (36)  =  2 . *  ( 1 . +POISSON) 

WRITE  (6,1909)  (AMAT (I) , 1=1, 36) 

1909  FORMAT  ('  COMPLIANCE  MATRIX  =' , 1P6E12 . 5/20X, 6E12 . 5/20X, 6E12 . 5/ 
1  2 OX, 6E12 . 5/20X, 6E12 .5/20X, 6E12 .5) 

END 


SUBROUTINE  CFCOMP(IBUG) 

C  CONSTRUCT  THE  COMPLIANCE  MATRIX 

C  1,1  =  1111;  1,2  =  1122;  1,3  =  1112;  .  .  .  3,3  =  1212 

INCLUDE  ' $DISK3 : [ SEAMAN . CRACK] KCKCOM . FOR' 

C 

DO  908  1=1,36 
908  AMAT(I)  =  0. 

AMAT(l)  =  1 .+PI/SAREA*SUM(1, 1) 

AMAT ( 2 )  =  -POISSON+PI/SAREA*SUM (1,2) 

AMATO)  =  -POISSON 

AMAT (4)  =  P I / SAREA*SUM ( 1 , 3) 

AMAT (7)  =  -POISSON+P I / S AREA* SUM (2,1) 

AMAT ( 8 )  =  1 . +PI/ S AREA* SUM (2,2) 

AMAT (9)  =  -POISSON 

AMAT (10)  =  PI/SAREA*SUM(2, 3) 

AMAT (13)  =  -POISSON 
AMAT (14)  =  -POISSON 
AMAT (15)  =  1. 

AMAT (19)  =  PI/SAREA*SUM(3, 1) 

AMAT (20)  =  PI/ S AREA*  SUM (3,2) 

AMAT (22 )  =  2 . * (1 ,+POISSON) +PI/SAREA*SUM(3, 3) 

AMAT (29)  =  2 .* (1 ,+POISSON) 

AMAT (36)  =  2 .* (1 ,+POISSON) 

WRITE  (6,1908) 

1908  FORMAT  ('  ************  COMPLIANCE  AND  STIFFNESS  MATRICES  FOR', 

1  '  INTERACTING  CRACKS,  ACCORDING  TO  KACHANOVS  METHOD  ********') 

WRITE  (6,1909)  (AMAT (I) , 1=1, 36) 

1909  FORMAT  ('  COMPLIANCE  MATRIX  =  ' , 1P6E12 . 5, '  FROM  CRCOMP  909' 

1  /20X, 6E12.5/20X, 6E12.5/20X, 6E12.5/20X, 6E12.5/20X, 6E12 .5) 

END 


C-2  7 


SUBROUTINE  CRFORC(IBUG) 

C  COMPUTE  THE  STRESSES  DUE  TO  THE  EXTERNAL  LOADING 
C 

INCLUDE  ' $DISK3 : [SEAMAN . CRACK] KCKCOM. FOR' 

C 

C  FILL  THE  FORCE  ARRAY  WITH  STRESSES  ON  EACH  C RACK 
C  FORCE (  +1)  IS  THE  NORMAL  COMPONENT,  FORCE (  +2)  THE  SHEAR  COMP. 

DO  200  NC  =  1 , NCRACK 

FORCE (2*NC-2+l)  =  SXX*COSTH (NC) **2  +  2 . *SXY*COSTH (NC) *SINTH (NC)  + 
1  SYY*SINTH (NC) **2 

FORCE (2*NC-2+2)  =  - (SXX-SYY) *SINTH (NC) *COSTH (NC)  + 

1  SXY* (COSTH (NC) **2-SINTH (NC) **2) 

200  CONTINUE 

KROW  =  2* NCRACK 
IF  ( IBUG  .EQ.  0)  GO  TO  210 
WRITE  (6,1504)  SXX, SYY, SXY 
WRITE  (6,1507)  (FORCE (I) , 1=1, KROW) 

1504  FORMAT  ('  FORC-205  FORCE  VECTOR:  NORMAL  AND  SHEAR  STRESS  ACTING  ' 
1  'ON  EACH  CRACK  FOR  INPUT  SXX,  SYY,  SXY=' , 1P3E10 . 3/10X, 5X, 

1  'NORMAL' , 10X, ' SHEAR' ) 

1507  FORMAT  ( 10X, 1P2E15 . 7 ) 

210  CONTINUE 
END 


no  oooooo 


SUBROUTINE  CRK1  (IBUG) 

DETERMINATION  OF  TRACTIONS  ON  THE  CRACK  FACES  IN  PREPARATION  FOR 
COMPUTING  THE  K1C  VALUES. 

INCLUDE  ' $DISK3 : [SEAMAN. CRACK [KCKCOM. FOR' 

CONSTRUCT  THE  ARRAYS  (SP  AND  ST)  OF  STRESSES  AT  POINTS  ALONG 
THE  NC  CRACK. 

DO  60C  NC  =  1 , NCRACK 
II  =  2*  (NC-1) 

DXEL  =  (XCRACK (NC,  3 ) -XCRACK (NC, 1 ) ) / 2 . 

DYEL  =  ( YCRACK (NC, 3) -YCRACK (NC, 1) ) /2 . 

DO  520  N  =  1 , NINTRV 
SUMN(N)  =  0. 

SUMT(N)  =  0. 

520  CONTINUE 

DO  580  J  =  1, NCRACK 
JJ  =  2* ( J-l) 

IF  (NC  .EQ.  J)  GO  TO  580 

SINA  =  SINTH (NC) *COSTH (J) -SINTH (J) *COSTH (NC) 

COSA  =  COSTH (NC) *COSTH ( J) +SINTH (NC) *SINTH ( J) 

SINA2  =  SINA* *2 
COSA2  =  COSA**2 
SINCOSA  =  SINA*COSA 

IF  (IBUG  .GT.  2)  WRITE  (6,1510)  NC, J, SINA, COSA 
1510  FORMAT  ('  Kl-510  NC=',I4,'  J=',I4,'  SINA, COSA=' , 1P2E12 . 5) 

DO  570  N  =  1, NINTRV 

XI  -  COS((2*(NINTRV+l-N)-l)*PI/(2.*NINTRV)) 

XID(N)  =  XI 

XL  =  XCR(NC)  +  XI*DXEL 
YL  =  YCR(NC)  +  XI*DYEL 

CALL  -TRACTN-  TO  DETERMINE  THE  TRACTIONS  ON  THE  CRACK  FACES 
CALL  TRACTN  (LS,  NC,  J,  XL,  YL,  XCR  (J)  ,  YCR  (J)  ,  SINTH,  COSTH,  ELCRAK  (J)  , 

1  SXXN, SYYN, SXYN, SXXT, SYYT, SXYT, IBUG) 

C 

IF  (IBUG  .GT.  2)  WRITE  (6,1568)  NC, J, N, SXXN, SYYN, SXYN, SXXT, SYYT, 
1  SXYT 

1568  FORMAT  ('  Kl-568  NC, J, N=' , 314, '  SXXN, Y, XY=' , 1P3E12 . 5, 

1  '  SXXT, Y, XY=' , 3E12 . 5) 

SNN (NC, J, N)  =  SXXN*COSA2  +  SYYN*SINA2  +  2 . *SXYN*SINCOSA 
STN (NC, J, N)  =  - (SXXN-SYYN) *SINCOSA  +  SXYN* (COSA2-SINA2) 

SNT (NC, J, N)  =  SXXT*COSA2  +  SYYT*SINA2  +  2 . *SXYT*SINCOSA 
STT (NC, J, N)  =  - (SXXT-SYYT) ‘SINCOSA  +  SXYT* (COSA2-SINA2) 

IF  (IBUG  .GT.  2)  WRITE  (6,1540)  SNN (NC, J, N) , STN (NC, J, N) , 

1  SNT  (NC,  J,  N)  ,  STT  (NC,  J,  N) 

1540  FORMAT  ('  Kl-540  SNN, STN=' , 1P2E12 . 5, '  SNT, STT=' , 2E12 . 5) 

SUMN(N)  =  SUMN(N)  +  SNN (NC, J, N) *PL ( 2  * J-l )  +  SNT (NC, J, N) *PL (2*J) 
SUMT(N)  =  SUMT(N)  +  STN (NC, J, N) *PL (2* J-l)  +  STT  NC, J, N) *PL (2 * J) 
570  CONTINUE 
580  CONTINUE 

DO  585  N  =  1, NINTRV 

SP (NC, N)  =  FORCE (2*NC-1)  +  SUMN (N) 

ST (NC, N)  =  FORCE (2*NC)  +  SUMT (N) 

585  CONTINUE 
590  CONTINUE 
6 00  CONTINUE 
END 


SUBROUTINE  CRK2 {NC, IBUG) 

C  COMPUTE  THE  CRACK  TIP  STRESS  INTENSITY  FACTORS 

C 

INCLUDE  ' $DISK3 : [SEAMAN .CRACK] KCKCOM. FOR' 

C 

WRITE  (6,1605)  NC 

1605  FORMAT  (/'  **********************', 40X, '**********************' / 

1  '  CRACK  TIP  STRESS  INTENSITIES  AND  CRACK  SHAPES  FOR 

2  'CRACK  NO.', 12,'  WRITTEN  BY  CRK2 ' / ) 

AK1  (NC,  1)  =  0  . 

AK1 (NC, 2)  =  0 . 

AK2  (NC,  1)  =  0. 

AK2 (NC, 2 )  =  0. 

DO  660  N  =  1 , NINTRV 

AK1  (NC,  1)  =  AK1  (NC,  1)  +  (1 . +XID  (N)  )  *SP  (NC,  N) 

AK1 (NC, 2)  »  AK1(NC,2)  +  ( 1 . -XID (N) ) *SP (NC, N) 

AK2  (NC,  1)  =  AK2  (NC,  1)  +  ( 1 . +XID  (N)  )  *ST  (NC,  N) 

AK2  (NC,  2 )  =  AK2  (NC,  2)  +  <1 . -XID  (N)  )  *ST  (NC,  N) 

IF  (IBUG  .GE.  2)  WRITE  (6,1567)  NC, N, AK1 (NC, 1) , AK1 (NC, 2) , 

1  AK2 (NC, 1) , AK2 (NC,2) , XID (N) , SP (NC,N) , ST(NC,N) 

1567  FORMAT  ('  K2-567  NC,N=',2I4,'  Kl=' , 1P2E10 . 3, '  K2=',2E10.3, 

1  /'  XID=' , E10 . 3, '  SP, ST=' , 2E10 . 3) 

660  CONTINUE 

SQPIL  =  SQRT (PI*ELCRAK (NC) ) 

AK1 (NC, 1 )  =  AK1 (NC, 1) /NINTRV*SQPIL 
AK1(NC,2)  =  AK1 (NC, 2) /NINTRV*SQPIL 
AK2 (NC, 1 )  =  AK2 (NC, 1) /NINTRV*SQPIL 
AK2 (NC, 2)  -  AK2 (NC, 2) /NINTRV* SQPIL 
PAVG  *  SXX*COSTH (NC) **2+SYY*SINTH (NC) **2 
1  +2 .  *SXY*SINTH (NC) *COSTH (NC) 

AK10  =  SQPIL*PAVG 

WRITE  (6,  1702)  AK1  (NC,  1)  ,  AK1  (NC,  2)  ,  AK10,  AK2  (NC,  1)  ,  AK2  (NC,  2 ) 

1702  FORMAT  ('  K  FACTORS' , 3X, ' RIGHT  TIP',3X,'  LEFT  TIP' , 4X, ' ISOLATED' , 

1  '  WRITTEN  BY  CRK2' / 

2  '  K1  =' , 1P3E12 .5/ 

3  '  K2  =' , 2E12 . 5) 

IF  (ABS(AKIO)  .LT.  l.E-10)  RETURN 
AK11  =  AK1 (NC, 1) /AK10 
AK12  =  AK1 (NC, 2) /AK10 
WRITE  (6,1704)  AK11 , AK12 
1704  FORMAT  ('  K1/K10' , 3H' S-, 1P2E12 . 5) 

END 


1 


SUBROUTINE  CRNINT  ( IBUG, SINALF, COSALF) 

C  COMPUTE  COMPLIANCE  MATRIX  FOR  NON-INTERACTING  CRACKS 

C 

INCLUDE  ' SDISK3 : [SEAMAN . CRACK] KCKCOM . FOR' 

C 

SUM { 1 , 1)  =  0 . 

SUM (2,1)  =  0. 

SUM (3, 1)  =  0. 

DO  10  NC  =  1 , NCRACK 

SUM (1,1)  =  SUM (1, 1) +ELCRAK (NC) **2* (COSTH (NC) *COSALF  + 

1  SINTH (NC) *SINALF) **2 

SUM (2,1)  =  SUM (2. 1) +ELCRAK (NC) * *2* (SINTH (NC) *COSALF  - 
1  COSTH (NC) * SINALF) **2 
SUM (3,  1)  =  SUM(3, 1) +ELCRAK(NC) **2 
10  CONTINUE 

DO  50  1=1,36 
50  AMAT(I)  =  0. 

AMAT(l)  =  1 .+PI/SAREA*SUM(1, 1) 

AMAT (2)  =  -POISSON 
AMAT ( 3 )  =  -POISSON 
AMAT ( 4 )  =  0 . 

AMAT (7)  =  -POISSON 

AMAT ( 8 )  =  1 .+PI/SAREA*SUM(2, 1) 

AMAT (9)  =  -POISSON 
AMAT (10)  =  0. 

AMAT (13)  =  -POISSON 
AMAT (14)  =  -POISSON 
AMAT (15)  =  1. 

AMAT (19)  =  0. 

AMAT (20)  =  0. 

AMAT (22)  =  2 . * (1 . +POISSON) +PI/SAREA*SUM (3, 1) 

AMAT (29)  =  2 . * ( 1 . +POISSON) 

AMAT (36)  =  2 . * (1 .+POISSON) 

WRITE  (6,1908) 

1908  FORMAT  ('  *****  COMPLIANCE  AND  STIFFNESS  MATRICES  FOR  NONINTERA' , 
1  'CTING  CRACKS  ****»') 

WRITE  (6,1909)  (AMAT (I) ,1=1,36) 

1909  FORMAT  ('  COMPLIANCE  MATRIX  =' , 1P6E12 . 5, '  FROM  CRNINT  909' 

1  /20X, 6E12 . 5/20X, 6E12 .5/20X, 6E12 .5/20X, 6E12.5/20X, 6E12.5) 

END 


ooo  oo  ono  oooo 


SUBROUTINE  CRSHAP  (NC, IBUG) 

COMPUTATION  OF  THE  CRACK  SHAPES 

INCLUDE  ' SDISK3 : [ SEAMAN. CRACK ] KCKCOM. FOR' 

COMPUTATION  OF  THE  ELLIPTIC  SHAPE  COEFFICIENTS 
SQPIL  =  SQRT (PI*ELCRAK (NC) ) 

SHAPEA (NC,  1 )  =  0 . 

SHAPEB (NC, 1)  =  0. 

SHAPEA (NC, 2)  =  0. 

SHAPEB (NC, 2)  =  0. 

SHAPES (NC,1)  =  2 . *PL (2*NC-1)  -  0 . 5* (AKI (NC, 1) +AK1 (NC, 2) ) /SQPIL 
IF  (ABS  (AKI  (NC,  1) +AK1  (NC,2)  )  .  LT  .  1.)  GO  TO  665 

SHAPEA (NC, 1)  =  (AKI (NC, 1 )  -  AKI (NC, 2 ) ) / ( 4 . *PL (2 *NC- 1 ) *SQPIL 
1  -  AKI (NC, 1)  -  AKI (NC, 2) ) 

SHAPEB (NC, 1 )  =  2 .* (AKI (NC, 1)  +AK1(NC,2)  -  2 . *PL (2*NC-1) *SQPIL)  / 
1  (4 . *PL (2*NC-1) *SQPIL  -  AKI (NC, 1)  -  AK1(NC,2)) 

665  SHAPES(NC,2)  =  2.*PL(2*NC)  -  0 . 5* (AK2 (NC, 1) +AK2 (NC, 2) ) /SQPIL 
IF  (ABS  (AK2  (NC,  1) +AK2  (NC,  2)  )  .LT.  I.)  GO  TO  670 
SHAPEA(NC,2)  =  (AK2 (NC, 1)  -  AK2 (NC, 2 ) ) / ( 4 . *PL (2 *NC) *SQPIL 
1  -  AK2  (NC,  1 )  -  AK2  (NC,  2)  ) 

SHAPEB (NC, 2)  =  2 . * (AK2 (NC, 1)  +  AK2(NC,2)  -  2 . *PL (2*NC) *SQPIL) / 

1  ( 4 . *PL (2*NC) *SQPIL  -  AK2(NC,1>  -  AK2(NC,2)) 


NORMAL  AND  SHEAR  OPENING  AREA 


670  CONTINUE 

WRITE  (6,1670)  (SHAPES (NC, I) , SHAPEA (NC, I) , SHAPEB (NC, I) , 1=1,2) 

1670  FORMAT  (/'  SHAP  670  SHAPE  FUNCTIONS  IN  THE  FORM  4L/E*S ( 1+Ax+Bx2 ' , 

1  ') Ellipse'/'  FACTORS  FOR  OPENING  ARE  S,  A,  B  =',1P3E11.3/ 

2  '  FACTORS  FOR  SHEARING  ARE  S,  A,  B  =',3E11.3) 

*******  RETURN  IN  CASE  PRINTING  IS  NOT  REQUESTED.  ******** 

IF  (IBUG  .LE.  0)  RETURN 

AREA  =  2 . *PI*ELCRAK (NC) * SHAPES (NC, 1) /EMOD* (1 . +0 . 25*SHAPEB (NC, 1) ) 
AREAT  =  2 . *PI*ELCRAK (NC) * SHAPES (NC, 2) /EMOD* (1 . +0 . 25*SHAPEB (NC, 2 ) ) 
DISPLACEMENTS  OF  CRACKS 

DXX  =  0.5*  (AREA*COSTH  (NC)  **2  +  AP.EAT*SINTH  (NC)  *COSTH  (NC)  ) 

DYY  -  0.5* (AREA*SINTH (NC) **2  -  AREAT*SINTH (NC) *COSTH (NC) ) 

DXY  =  0 . 5*AREA*SINTH (NC) *COSTH (NC)  -  0 . 25*AREAT* (COSTH (NC) **2 
1  -SINTH (NC) **2) 

WRITE  (6,1677)  NC,  AREA,  AREAT,  DXX,  DYY,  DXY 
1677  FORMAT  (/'  SHAP  677  OPENINGS  FOR  CRACK  ',13,'  AREA, AREAT=' , 

1  1P2E11 . 3, '  DXX,  DYY,  DXY=' , 3E11 . 3) 


COMPUTATION  OF  CRACK  SHAPES 


NPLOT  =  20 

OPISO  =  4 . *ELCRAK (NC) /EMOD*PAVG 
VSHEAR  -  0 . 

XOL  =  -1. 

VNORM  =*  0. 

WRITE  (6,1688)  NC 

1688  FORMAT  ('  SHAP  688  CRACK  SHAPE  FOR  CRACK  NO  ' , I3/4X,  '  N' ,  8X,  '  X/L' , 
1  9X, ' BN' , 9X, ' BT' , 6X, 'BNISO' , 5X, 'VSHEAR' , 6X, 'VNORM' , 4X, ' ELLIPSE' ) 

N  =  1 
BNISO  =  0. 

WRITE  (7,1703)  N, XOL, VSHEAR, VNORM, BNISO 
WRITE  (6,1703)  N, XOL, VSHEAR, VNORM, BNISO, VSHEAR, 

1  VNORM, BNISO 


DO  680  N  =  2 , NPLOT 

RATIO  =  FLOAT (M-1) /FLOAT (NPLOT) 

XOL  =  COS ( (1. -RATIO) *PI) 

ELLIP  =  SQRT ( 1 . -XOL*  *2 ) 

VNORM  =  (1.+  SHAPEA(NC, 1) *XOL  +  SHAPEB (NC, 1 ) *XOL*  *2 ) *ELLIP 
VSHEAR  =  (1.+  SHAPEA (NC, 2 ) *XOL  +  SHAPEB (NC, 2 ) *XOL* *2 ) *ELLIP 
BN  =  4 . *ELCRAK (NC) /EMOD*SHAPES (NC, 1 ) *VNORM 
BT  =  4 . *ELCRAK (NC) /EMOD*SHAPES (NC, 2 ) *VSHEAR 
BNISO  =  OP ISO*ELLIP 

WRITE  (6,1703)  N, XOL, BN, BT, BNISO, VSHEAR, VNORM, 

1  ELLIP 

1703  FORMAT  ( 15 , 1P10E11 . 3 ) 

WRITE  (7,1703)  N, XOL, VSHEAR, VNORM, ELLIP 
680  CONTINUE 

N  =  NPLOT+1 
XOL  =  1. 

VSHEAR  =  0 . 

VNORM  =  0 . 

WRITE  (6,1703)  N, XOL, VSHEAR, VNORM, VNORM, VSHEAR, 

1  VNORM, VNORM 

WRITE  (7,  1703)  N,  XOL,  VSHEAR, VNORM, VNORM 
WRITE  (6,1684) 

1684  FORMAT  ('  DEFINITIONS:  BN  and  BT  are  normal  (+  in  opening)  and 

1  shear  (+  CCW)  displacement  for  the  cracks  (cm) ,  (FROM  SHAP  684) ' 

2  /14X, 'BNISO  is  the  normal  opening  for  ', 

3  'an  isolated  crack  (cm) ' /14X, ' VNORM  and  VSHEAR  are  normalized', 

4  '  opening  and  shearing  distortion' /14X, ' ELLIPSE  is  the  ' 

5  'normalized  elliptic  distortion') 

END 


SUBROUTINE  CRSTIF  ( I BUG, S INALF, COS ALF) 

C  Computation  of  SUM  Lk~2  (NkTk  +  TkNk)  for  each  unit  stress  loading. 

C 

INCLUDE  ' $DISK3 : ( SEAMAN. CRACK] KCKCOM. FOR' 

C 

C  Construct  STRESS  TENSORS  in  X-Y  plane  corresponding  to  Unit 

C  Stresses  Sll,  S22,  and  S12  in  the  orthotropic  axes. 

C  CCW  Rotation  of  the  STRESS  back  to  X-Y  axes  by  an  angle  -ALF. 

DO  900  I LOAD  =1,3 
IF  ( I LOAD  . GT .  1)  GO  TO  800 
C  TENSOR  FOR  THE  Sll  LOADING 

SXX  =  COSALF**2 
SYY  =  SINALF**2 
SXY  =  S INALF *COS ALF 
GO  TO  820 

800  IF  (I LOAD  .GT.  2)  GO  TO  810 
C  TENSOR  FOR  THE  S22  LOADING 

SXX  =  S INALF* *2 
SYY  =  COSALF**2 
SXY  =  -S INALF *COS ALF 
GO  TO  820 

C  TENSOR  FOR  THE  S12  LOADING 

810  CONTINUE 

SXX  =  -2 . *S INALF *COS ALF 
SYY  =  2 . *S INALF *COSALF 

SXY  =  COSALF**2-SINALF*  *2 
820  CONTINUE 

C - 

C  FILL  THE  FORCE  ARRAY  WITH  STRESSES  ON  EACH  CRACK 
C  FORCE (  +1)  IS  THE  NORMAL  COMPONENT,  FORCE (  +2)  THE  SHEAR  COMP. 

C - 

CALL  CRFORC  (IBUG) 

C 

DO  840  I  -  1 , KROW 
840  PL (I)  =  FORCE (I) 

KEND  =  4*NCRACK**2 
DO  845  K  =  1, KEND 
845  BMAT(K)  =  AMAT(K) 

C - 

C  COMPUTE  THE  STRESSES  ON  EACH  MICROCRACK,  ACCOUNTING  FOR  BOTH  THE 

C  APPLIED  LOAD  AND  ALL  OTHER  CRACKS.  -FORCE-  IS  APPLIED,  -PL-  IS 

C  THE  STRESS  (NORMAL,  SHEAR)  ON  EACH  CRACK. 

C - 

CALL  SIMQ (BMAT, PL, 2 *NCRACK, KSTOP ) 

IF  (IBUG  .GT.  0)  THEN 

WRITE  (6,1841)  I LOAD,  SXX,  SYY,  SXY 
1841  FORMAT  (/'  STIF  841  I LOAD  =',I3,'  SXX, SYY, SXY=' , 1P3E15 . 8 ) 

WRITE  (6,1845)  (FORCE (I) , 1=1, KROW) 

1845  FORMAT  ('  STIF  845  FORCE  (N,  TAU)  =  ',1P7E15.8) 

WRITE  (6,1843)  (PL(I) ,1=1, KROW) 

1843  FORMAT  ('  STIF  843  P  (N,  TAU)  =  ',1P7E15.8) 

WRITE  (6,1848) 

1848  FORMAT  ('  (STIF)  AMAT  MATRIX'  ) 

DO  850  I  =  1 , KROW 
KK  =  ( I  —  1 )  *KROW 

WRITE  (6,1850)  (AMAT (KK+L) ,L=1, KROW) 

850  CONTINUE 

1850  FORMAT  (1P7E15 . 8/ (5X, 7E15 . 8) ) 

ENDIF 
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T1  IS  STRESS  IN  THE  1  DIRECTION,  T2  IN  THE  2  DIRECTION 
IN  THE  ORTHOTRCP IC  COORDINATES 


SUM ( 1 , I LOAD )  =  0. 

SUM (2 , I LOAD )  =  0 . 

SUM (3, I LOAD)  =  0. 

DO  860  NC  =  1 , NCRACK 

COSTHA  =  COSTH (NC) *COSALF  +  SINTH (NC) *SINALF 
SINTHA  =  SINTH (NC) *COSALF  -  COSTH (NC) *SINALF 
T1  =  PL (2*NC-2+l) *COSTHA  -  PL (2*NC-2+2) *SINTHA 
T2  =  PL (2*NC-2+l) *SINTHA  +  PL(2*NC-2+2) *COSTHA 
EL2  =  ELCRAK(NC) **2 

IF  ( I BUG  .GT.  0)  WRITE  (6,1852)  NC,  T1,T2 
1852  FORMAT  ('  STIF  852  CRACK  NO', 12,'  STRESSES  IN  1  AND  2  DIRECTIONS,' 
1  '  T1  AND  T2  =' , 1P2E12 . 5 ) 


SUM  THE  STRESS  QUANTITIES  FOR  ALL  CRACKS  TO  CONSTRUCT  THE 
SUM  La2  (NkTk  +  TkNk)  terms 


DSUM1  =  EL2*Tl*COSTHA 
SUM ( 1 , ILOAD)  =  SUM ( 1, ILOAD)  +  DSUM1 
DSUM2  =  EL2*T2* SINTHA 
SUM(2, ILOAD)  =  SUM(2, ILOAD)  +  DSUM2 
DSUM3  =  EL2* (T2*COSTHA  +  T1*SINTHA) 

SUM(3,  ILOAD)  =  SUMO,  ILOAD)  +  DSUM3 

IF  ( IBUG  .GT.  1)  WRITE  (6,1858)  ILOAD, DSUM1, DSUM2, DSUM3, 

1  (SUM (I, ILOAD) ,1=1,3) 

1858  FORMAT  ('  STIF  858  IL0AD=',I3,'  DSUMs  =' , 1P3E12 . 5, '  SUM=' , 
1  1P3E12.5) 

860  CONTINUE 
900  CONTINUE 
END 
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SUBROUTINE  CRSTRS(IBUG) 

COMPUTATION  OF  THE  STRESSES  ON  EACH  MICROCRACK 
INCLUDE  ' SDISK3 : [SEAMAN .CRACK] KCKCOM . FOR' 


BEGIN  LOOP  TO  COMPUTE  THE  L(I,J)  MATRIX 


I  IS  THE  INDICATOR  FOR  THE  MICROCRACK  OF  INTEREST 
J  IS  ANOTHER  MICROCRACK 
NCR2  =  2  *NCRACK 
ANINTR3  =  3*NINTRV 
DO  500  NC  -  1 ,  NCRACK 
IN  =  2*NC-1 
IT  =  2*NC 

DX  =  (XCRACK (NC, 3) -XCRACK (NC, 1) } /NINTRV 
DY  =  (YCRACK(NC, 3) -YCRACK(NC, 1) ) /NINTRV 
DO  450  J  *  1, NCRACK 
JN  =  2* J-l 
JT  =  2*J 

IF  (NC  .NE.  J)  GO  TO  220 
AMAT ( (JN-1) *NCR2+IN)  =  1. 

AMAT ( (JT-1) *NCR2+IN)  =  0. 

AMAT ( (JN-1) *NCR2+IT)  -  0. 

AMAT ( (JT-1) *NCR2+IT)  =  1. 

GO  TO  450 
220  CONTINUE 

SINA  =  SINTH(NC) *COSTH ( J) -SINTH ( J) *COSTH (NC) 

COSA  =  COSTH(NC) *COSTH ( J) +SINTH (NC) *SINTH(J) 

SINA2  =  SINA* *2 
COSA2  =  COSA* *2 
SINCOSA  =  SINA*COSA 
SUMXXN  =  0. 

SUMYYN  -  0. 

SUMXYN  =  0. 

SUMXXT  *  0 . 

SUMYYT  =  0. 

SUMXYT  =  0. 

NINTRV1  =  NINTRV+1 
DO  300  N  -  1 , NINTRV1 
XL  =  XCRACK  (NC,  1)  +  (N-l)  *DX 
YL  =  YCRACK (NC, 1) + (N-l) *DY 

CALL  -TRACTN-  TO  DETERMINE  THE  TRACTIONS  ON  THE  CRACK  FACES 
CALL  TRACTN  (LS, NC, J, XL, YL, XCR (J) , YCR (J) , SINTH, COSTH, ELCRAK (J) , 

1  SXXN, SYYN, SXYN, SXXT, SYYT, SXYT, IBUG) 

C 

IF  (IBUG  .GT.  2)  WRITE  (6,1297)  NC, J, N, SXXN, SYYN, SXYN, SXXT, SYYT, 
1  SXYT 

1297  FORMAT  ('  STRS  297  NC, J,N=' , 314, '  SXXN, Y, XY=' , 1P3E12 . 5 , 

1  '  SXXT, Y, XY=' , 3E12 . 5) 

COEF  -  1. 

IF  (N  .GT.  1  .AND.  N  . LT .  NINTRV1)  COEF  -  2 . +2 . *MOD (N-l , 2 ) 

SUMXXN  -  SUMXXN  +  COEF* SXXN 
SUMYYN  =  SUMYYN  +  COEF* SYYN 
SUMXYN  =  SUMXYN  +  COEF*SXYN 
SUMXXT  -  SUMXXT  +  COEF* SXXT 
SUMYYT  -  SUMYYT  +  COEF* SYYT 
SUMXYT  -  SUMXYT  +  COEF*SXYT 


300  CONTINUE 

C  CONSTRUCT  THE  STRESSES  ON  THE  NC  CRACK  IN  THE  ORIENTATION  OF  THE 
C  NC  CRACK:  STRNN  AND  STRNT  -  NORMAL  STRESS  ON  NC  FROM  NORMAL  AND 
C  SHEAR  ON  J,  AND  STRTN  AND  STRTT  -  SHEAR  ON  NC  FROM  NORMAL  AND 

C  SHEAR  ON  J. 

STRNN  =  ( SUMXXN*COSA2  +  SUMYYN*SINA2  +  2 . *SUMXYN*SINCOSA) / ANINTR3 
STRTN  =  (- (SUMXXN-SUMYYN) *SINCOSA  +  SUMXYN* (COSA2 -S INA2 ) ) / ANINTR3 
STRNT  =  ( SUMXXT*COSA2  +  SUMYYT*SINA2  +  2 . *SUMXYT*SINCOSA) / ANINTR3 
STRTT  =  {- (SUMXXT-SUMYYT) *SINCOSA  +  SUMXYT* (COSA2-SINA2 ) ) /ANINTR3 
C  CONSTRUCT  THE  -A-  MATRIX:  A (EFFECT, CAUSE) 

AMAT ( (JN-1) *NCR2+IN)  =  -STRNN 
AMAT ( (JT-1) *NCR2+IN)  =  -STRNT 
AMAT ( (JN-1) *NCR2+IT)  =  -STRTN 
AMAT ( (JT-1) *NCR2+IT)  =  -STRTT 
450  CONTINUE 
500  CONTINUE 

KROW  =  2  *NCRACK 

KEND  =  4  *NCRACK*  *2 

IF  ( IBUG  .EQ.  0)  GO  TO  508 

WRITE  (6,1502) 

DO  505  I  =  1 , KROW 
KK  =  ( 1-1 ) *KROW 

WRITE  (6,1503)  (AMAT (KK+L) , L=l, KROW) 

505  CONTINUE 
WRITE  (6,1500) 

DO  506  I  =  1 , KROW 

WRITE  (6,1503)  (TL (I, J) , J=l, KROW) 

506  CONTINUE 

508  CONTINUE 

DO  509  K  =  1 , KEND 

509  BMAT(K)  =  AMAT (K) 

1500  FORMAT  (/'  STRS  500,  -TL-  MATRIX') 

1502  FORMAT  (/'  STRS  502,  -AMAT-  MATRIX') 

1503  FORMAT  (1P7E15 . 8/ (5X,7E15.8) ) 

C - 

C  SOLUTION  FOR  STRESSES  ON  EACH  MICROCRACK,  CONSIDERING  BOTH 
C  THE  APPLIED  STRESS  FIELD  AND  ALL  OTHER  CRACKS 

C  PL ( 1 )  =  NORMAL  STRESS  AND  PL (2)  =  SHEAR  STRESS 

C - 

DO  510  I  =  1 , KROW 

510  PL ( I )  =  FORCE (I) 

CALL  SIMQ (BMAT, PL, 2*NCRACK,  KSTOP) 

IF  (IBUG  .GT.  0)  WRITE  (6,1550)  (PL (I) , 1=1, KROW) 

1550  FORMAT  ('  STRS  550,  AVERAGE  FORCE  ACTING  ON  EACH  CRACK:  ', 

1  'P  (N,  TAU) ' /19X, 'NORMAL' , 10X, ' SHEAR' / (10X, 1P2E15 . 8)  ) 

END 
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SUBROUTINE  MINV ( A, N, D, L, M) 

IMPLICIT  REAL*8 (A-H,0-Z) 

DIMENSION  A ( 1 ) ,L(1) ,M(1) 

SEARCH  FOR  LARGEST  ELEMENT 
D=l. 

NK=-N 

DO  80  K=1 , N 
NK=NK+N 
L (K) =K 
M  (K)  =K 
KK=NK+K 
BIGA=A (KK) 

DO  20  J=K,N 
IZ=N* (J-l) 

DO  20  I=K,N 
IJ=IZ+I 

IF  (ABS (BIGA) -ABS (A(IJ) ) )  15,20,20 
BIGA=A ( I J) 

L(K)=I 
M (K) =J 
CONTINUE 

INTERCHANGE  ROWS 
J=L (K) 

IF  (J-K)  35,35,25 

KI=K-N 

DO  30  1=1, N 

KI=KI+N 

HOLD=-A(KI) 

JI=KI-K+J 

A(KI)=A(JI) 

A ( JI ) =HOLD 

INTERCHANGE  COLUMNS 
I=M(K) 

IF  (I-K)  45,45,38 
JP=N* (1-1) 

DO  40  J=1,N 
JK=NK+J 
JI=JP+J 
HOLD=-A(JK) 

A(JK)=A(JI) 

A ( JI ) =HOLD 

DIVIDE  COLUMN  BY  MINUS  PIVOT  (VALUE  OF  PIVOT  ELEMENT  IS 
CONTAINED  IN  BIGA) 

IF  (ABS (BIGA) -l.E-20)  46,46,48 
D=0  . 

WRITE  (6,446)  BIGA 

RETURN 

DO  55  1=1, N 

IF  (I-K)  50,55,50 

IK=NK+I 

A  ( IK)  =A  ( IK)  /  (-BIGA) 

CONTINUE 

REDUCE  MATRIX 
DO  65  1=1, N 
IK=NK+I 
HOLD=A ( IK) 

I  J=I-N 
DO  65  J=1,N 
I J=I J+N 


IF  (I-K)  60,65,60 

IF  (J-K)  62,65,62 
KJ=IJ-I+K 

A  ( I J) =HOLD*A(KJ) +A ( I J) 

CONTINUE 

DIVIDE  ROW  BY  PIVOT 
KJ=K-N 
DO  75  J=1 , N 
KJ=KJ+N 

IF  (J-K)  70,75,70 
A (KJ) =A(KJ) /BIGA 
CONTINUE 

PRODUCT  OF  PIVOTS  AND  REPLACE  PIVOT  BY  RECIPROCAL 
D=D*BIGA 
A (KK) =1 . /BIGA 
CONTINUE 

FINAL  ROW  AND  COLUMN  INTERCHANGE 

K=N 
K=K-1 

IF  (K)  150,150,105 
I=L(K) 

IF  (I-K)  120,120,108 
JQ=N* (K-l ) 

JR=N* (1-1) 

DO  110  J=1,N 
JK=JQ+J 
HOLD=A(JK) 

JI= JR+ J 
A ( JK) =-A(JI) 

A ( JI ) =HOLD 
J=M(K) 

IF  (J-K)  100,100,125 
KI=K-N 

DO  130  1=1, N 
KI=KI+N 
HOLD=A(KI) 

JI=KI-K+J 
A (KI ) =-A ( JI ) 

A ( JI ) =HOLD 
GO  TO  100 
RETURN 

FORMAT (/'  MINV  -  MATRIX  IS  SINGULAR,  BIGA  =  '1PE12.4/) 
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SUBROUTINE  S IMQ ( A, B, N, KS ) 

IMPLICIT  REAL*8  (A-H, O-Z) 

DIMENSION  A(l) , B (1) 

FORWARD  SOLUTION 
TOL=0 . 0 
KS=0 
JJ=-N 

DO  65  J=1 ,  N 
JY=J+1 
JJ=JJ+N+1 
BIGA=0 . 

IT=JJ-J 
DO  30  I=J,N 

SEARCH  FOR  MAXIMUM  COEFFICIENT  IN  COLUMN 
I J=IT+I 

IF  (ABS (BIGA) -ABS (A(IJ) ) )  20,30,30 
20  BIGA=A ( I J) 

IMAX=I 
30  CONTINUE 

TEST  FOR  PIVOT  LESS  THAN  TOLERANCE  (SINGULAR  MATRIX) 
IF  (ABS (BIGA) -TOL)  35,35,40 
35  KS=1 
RETURN 

INTERCHANGE  ROWS  IF  NECESSARY 
40  Il=J+N*(J-2) 

IT=IMAX-J 
DO  50  K=J,N 
I1=I1+N 
I2=I1+IT 
SAVE=A ( I 1 ) 

A (II) =A(I2) 

A ( 12 ) =SAVE 

DIVIDE  EQUATION  BY  LEADING  COEFFICIENT 
50  A(I1) =A(I1) /BIGA 
SAVE=B (IMAX) 

B (IMAX) =B (J) 

B (J) =SAVE/BIGA 

ELIMINATE  NEXT  VARIABLE 
IF  (J-N)  55,70,55 
55  IQS=N* ( J-l ) 

DO  65  IX=JY, N 

IXJ=IQS+IX 

IT=J-IX 

DO  60  JX=JY,N 

IXJX=N* (JX-1) +IX 

JJX=IXJX+IT 

6  0  A(IXJX) -A(IXJX) - (A(IXJ)  *A( JJX)  ) 

65  B(IX)«B<IX)-(B(J) *A(IXJ) ) 

BACK  SOLUTION 
70  NY=N-1 
IT=N*N 

DO  80  J=1 , NY 
IA=IT- J 
IB=N-J 
IC=N 

DO  80  K*1,J 

B  (IB)  =B  (IB  /  -A  (IA)  *B  ( IC) 

IA=*IA-N 
80  IC-IC-1 

C-40 


SUBROUTINE  TRACTN  (LS, NC, J, XL, YL, XCR, YCR, SINTH, COSTH, EL, 

1  SXXN, SYYN, SXYN, SXXT, SYYT, SXYT,  IBUG) 

C  ROUTINE  TO  COMPUTE  THE  TRACTIONS  ALONG  THE  LENGTH  OF  THE 

C  -NC-  CRACK  IN  RESPONSE  TO  THE  STRESSES  FROM  THE  -J-  CRACKS. 

C  CALLED  BY  KCRACK,  WRITTEN  BY  L.  SEAMAN  ON  JULY  29,  1987. 

IMPLICIT  REAL*8 (A-H,0-Z) 

DIMENSION  SINTH (10) , COSTH (10) 

X  =  (XL-XCR) *SINTH (J) - (YL-YCR) *COSTH (J) 

Y  =  (XL-XCR) *COSTH (J) + (YL-YCR) *SINTH (J) 

AL  =  SQRT ( (X-EL) **2+Y**2) 

GL  =  SQRT ( (X+EL) **2+Y**2) 

DL  =  SQRT (2 . * (X**2+Y**2-EL**2  + 

1  SQRT ( (X**2+Y**2+EL**2) **2  -  4 . * (X*EL) **2) ) ) 

IF  (IBUG  .GT.  2)  WRITE  (6,1010)  LS, NC, J, XL, YL, XCR, YCR, X, Y, 

1  SINTH (J) , COSTH (J) , SINTH (NC) , COSTH (NC) ,  SINA,  COSA,  EL, 

2  AL, DL, GL 

1010  FORMAT  ('  TRACTN  LS, NC, J=' , 313, '  XL, YL=' , 1P2E12 . 5, '  XCR,YCR=', 

1  2E12.5,'  X, Y=' , 2E12 .5/ '  SINJ, COSJ=' , 2E10 . 3, '  SINNC, COSNC=' , 

2  2E10.3,'  SINA, COSA=' ,2E10.3,'  EL=',E10.3/'  AL, DL, GL=' , 3E10 . 3 ) 

C  COMPUTE  Gi  =  Ii,  FACTORS  IN  THE  "STANDARD"  STRESS  FIELDS 

G1  =  4 . *EL**3* (GL-AL) / (DL* (AL+GL+DL) **2) 

G2  =  4 . *EL**2/ (DL* (AL+GL+DL) ) 

G3  =  2 . *EL**3* (GL-AL) / (AL*GL*DL**3) 

G4  =  2 . *EL**2* (AL+GL) / (AL*GL*DL**3) 

G5  =  0 . 5*EL**3* (3 . *AL*GL* (AL+GL) **2* (GL-AL) +DL**2* (GL**3-AL* *3 ) ) / 
1  ( (AL*GL) **3*DL**5) 

G6  »  0 .5*EL**2* ( (AL**3+GL**3) *DL**2  +  3 . *AL*GL* (AL+GL) **3) / 

1  ( (AL*GL) **3*DL**5) 

IF  (IBUG  .GT.  2)  WRITE  (6,1020)  G1,G2,G3,G4,G5,G6 
1020  FORMAT  ('  TRACTN  20  Gl, 2, 3*' , 1P3E22 . 5, '  G4, 5, 6=' , 3E12 . 5) 

C  COMPUTE  SXX,  SYY,  SXY,  THE  STRESSES  AT  THE  POINTS  X,Y  ON  THE  NC 
C  CRACK  CAUSED  BY  CONSTANT  NORMAL  AND  SHEARING  STRESSES  OF  UNIT 
C  INTENSITY  ON  THE  J  CRACK.  ORIENTATION  IS  WITH  RESPECT  TO  THE 
C  J  CRACK. 

SXXN  =  G2  -  8 . *Y**2*G4  +  8.*Y**4*G6 
SYYN  =  G2  +  4 . *Y**2*G4  -  8.*Y**4*G6 

SXYN  =  2.*(-Y*G3  +  X*Y*G4  +  4.*Y**3*G5  -  4 . *X*Y**3*G6) 

SXXT  =  2 . * (3 . *Y*G3  -  3 . *X*Y*G4  -  4.*Y**3*G5  +  4 . *X*Y**3*G6) 

SYYT  =  2 . * (-Y*G3  +  X*Y*G4  +  4.*Y**3*G5  -  4 . *X*Y**3*G6) 

SXYT  =  G2  -  8 . *Y**2*G4  +  8.*Y**4*G6 

IF  (IBUG  .GT.  2)  WRITE  (6,1030)  SXXN, SYYN, SXYN, SXXT, SYYT, SXYT 
1030  FORMAT  ('  TRACTN  30  SXXN, Y, XY=' , 1P3E12 . 5, '  SXXT, Y, XY=' , 3E12 . 5) 
RETURN 
END 
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